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ABSTRACT 


Future  spacecraft  designs,  and  in  particular  military  spacecraft,  may  incorporate  the  use  of 
synergetic  orbital  plane  change  maneuvers.  The  analysis  of  these  maneuvers  and  their 
optimality  is  an  area  in  which  much  work  has  been  done  but  only  a  few  questions  have 
been  answered.  This  thesis  discusses  the  theoretical  background  for  solving  the  optimal 
control  problem.  A  framework  is  set  forth  for  the  formulation  of  the  overall  problem 
which  must  be  solved.  Pontryagin’s  Maximum  Principle  is  applied  to  obtain  the  necessary 
conditions  for  maximizing  the  inclination  change  for  a  given  amount  of  propellant.  Effects 
of  a  heating  rate  constraint  imposed  by  the  thermal  protection  system  are  considered.  The 
Program  to  Optimize  Simulated  Trajectories  (POST)  is  used  to  obtain  results  for  the 
Maneuverable  Reentry  Research  Vehicle  (MRRV)  to  illustrate  certain  points.  Two 
characterizations  of  the  atmospheric  pass  are  analyzed  and  compared  to  previous  work, 
namely  Aerobang  and  Aerocruise.  A  discussion  on  the  limited  use  of  POST  as  a  direct 
method  of  analysis  is  also  included. 
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I.  INTRODUCTION 


Aeroassisted  orbital  plane  changes  offer  the  potential  of  significant  fuel  savings 
over  purely  exoatmospheric  propulsive  maneuvers.  Some  of  the  possible  applications  for 
this  maneuver  are: 

(1)  an  evasive  maneuver  for  a  satellite  being  pursued  by  an  enemy  AS  AT, 

(2)  an  ASAT  constellation  containing  many  ASAT’s  parked  in  some  benign 
looking  orbit  that  maneuver  when  needed  into  the  orbit  of  a  target  satellite, 

(3)  for  an  intelligence  satellite  to  gather  intelligence  at  a  veiy  low  altitude  or  to 
visit  a  specified  location  faster  than  an  exoatmospheric  satellite, 

(4)  for  a  reusable  ferry  between  a  repair  station  and  a  damaged  satellite  at 
different  inclinations, 

(5)  for  a  repair  satellite  to  perform  repairs  local  to  a  damaged  satellite  while 
being  stationed  at  a  repair  station. 

Some  of  these  orbital  transfer  applications  may  require  large  plane  changes  to  be  made  in 
a  single  pass.  Much  work  has  been  done  in  the  area  of  analyzing  possible  trajectories  for 
this  type  of  maneuver.  Three  distinct  trajectories  are  discussed  in  the  literature  as 
potential  optimal  solutions.  These  are  Aeroglide,  Aerocruise,  and  Aerobang.  All  of 
these  maneuvers  involve  a  deorbit  impulse,  an  aerodynamic  turn  in  the  atmosphere,  a 
reorbit  impulse  and  a  circularization  impulse.  The  difference  between  the  maneuvers  is 
in  the  aerodynamic  turn.  [Refs.  1-6] 

For  Aeroglide,  the  vehicle  performs  the  plane  change  solely  with  the  lift 
generated  by  the  vehicle  surfaces.  A  portion  of  the  lift  generated  is  normal  to  the  orbit 
plane  and  the  angle  of  attack  is  usually  close  to  (L/D)max.  In  order  to  generate  the 
necessary  lift,  the  vehicle  is  flown  at  high  velocities  and  low  altitudes  where  the 
atmosphere  is  dense  resulting  in,  sometimes,  excessively  high  heating  rates.  Aeroglide  is 
an  unconstrained  maneuver  in  terms  of  heating  rate,  meaning  that  no  heating  rate  limit  is 
taken  into  account.  The  only  fuel  expenditures  are  in  the  deorbit,  reorbit  and 
circularization  impulses  which  has  been  shown  to  be  significantly  less  than  the  purely 
propulsive  exoatmospheric  case  [Ref.  7],  The  tradeoff  between  Thermal  Protection 
System  (IPS)  mass  and  fuel  consumption  is  a  design  area  which  restricts  the  amount  of 
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pi  cine  change  which  can  be  achieved  by  an  Aeroglide  maneuver.  When  a  maneuver  is 
required  to  achieve  a  larger  plane  change  than  is  possible  with  the  unconstrained 
Aeroglide,  propulsive  power  must  be  used  to  control  the  heat  rate  during  the  turn. 
Aerocruise  and  Aerobang  are  two  schemes  for  doing  just  that. 

The  aerodynamic  turn  in  both  Aerocruise  and  Aerobang  consists  of  a  zero  thrust 
arc  from  atmosphere  entry  until  the  maximum  heating  rate  is  reached  and  a  propulsive 
arc  at  constant  heating  rate  followed  by  a  zero  thrust  arc  to  atmosphere  exit.  The 

difference  between  the  two  maneuvers  is  in  the  characteristics  of  the  heat  rate  control 
arc. 

Aerocruise  uses  a  throttled  thrust  to  cancel  drag  and  a  component  of  lift  to 
maintain  a  constant  altitude.  This  results  in  a  circular  arc  at  constant  velocity  and 
altitude  maintaining  a  constant  heating  rate.  Aerobang,  on  the  other  hand,  uses  zero 
thrust  or  maximum  thrust  and  modulates  the  angle  of  attack  and  angle  bank  to  vary 
altitude  (and  thus  density)  and  speed  while  maintaining  a  constant  heating  rate. 

The  question  then  arises,  “What  is  the  optimal  maneuver  to  maximize  the  orbital 
plane  change  of  a  circular  orbit  for  a  given  amount  of  propellant?”  The  answer  to  this 
question  remains  unknown.  This  thesis  will  build  directly  on  the  work  of  two  previous 
authors.  The  first  is  Thomas  P.  Spriesterbach  who  performed  a  comparison  of  two 
vehicles  designed  for  synergetic  maneuvers  and  their  performance  in  all  three  of  the 
above  profiles  [Ref.  2],  The  two  vehicles,  one  of  which  will  be  discussed  in  more  detail 
later,  are  the  Entry  Research  Vehicle  (ERV)  and  the  Maneuverable  Reentry  Research 
Vehicle  (MRRV).  In  his  analysis,  which  considered  only  at  the  atmospheric  portion  of 
the  maneuver,  he  concluded  that  the  MRRV  was  superior  to  the  ERV  in  performance 
and  that  the  Aerobang  maneuver  was  superior  to  Aerocruise  in  terms  of  plane  change  in 
the  atmospheric  portion  of  the  maneuver.  His  analysis  did  not,  however,  consider  the 
overall  maneuver  from  orbit  to  orbit.  John  C.  Nicholson,  who  used  the  Program  to 
Optimize  Simulated  Trajectories  (POST)  to  compare  the  three  profiles  for  the  orbit  to 
orbit  maneuver  using  the  ERV  [Ref.  1],  concluded  that  the  Aerobang  maneuver  was 
superior  in  inclination  change  to  Aerocruise.  Neither  of  these  works  included  theoretical 
justification  for  the  choice  of  the  maneuvers  or  an  overall  orbit  analysis. 
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This  thesis  extends  the  above  works  with  a  major  emphasis  placed  on  the 
theoretical  setup  of  the  optimal  control  problem.  It  is  important  to  note  form  the  start 
that  this  work  was  done  in  an  attempt  to  learn  about  the  optimal  trajectory  for  the 
aeroassisted  orbital  plane  change  from  a  low  earth  circular  orbit  to  another  low  earth 
circular  orbit  of  the  same  radius.  As  will  be  seen,  even  this  apparently  narrow  problem 
can  be  very  complex.  Although  later  there  will  be  numerical  analysis  of  two  particular 
feasible  trajectories  no  claim  is  being  made  here  as  to  the  optimality  of  any  particular 
trajectory.  However,  there  is  theoretical  data  available  which  narrows  the  search  when 
looking  for  the  optimal  trajectory.  It  is  difficult  to  solve  the  resulting  two  point  boundary 
value  problem  (TPBVP)  due  to  the  typical  small  radius  of  convergence'  however,  much 
can  be  learned  from  the  theory  which  can  aid  in  the  direct  methods  to  numerically  solve 
this  problem.  A  discussion  is  also  included  on  how  to  the  extend  this  thesis  to  obtain  the 
needed  knowledge  of  the  costate  values  and  allow  a  solution  to  the  TPBVP.  This  thesis 
develops  the  set  of  necessary  conditions  for  the  optimal  trajectory  and  discusses  a 
method  for  handling  the  small  radius  of  convergence. 
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II.  PROBLEM  STATEMENT  AND  MODELING 


For  the  purpose  of  this  thesis  the  term  optimal  will  be  used  in  the  sense  of  the 
largest  inclination  change  for  a  given  amount  of  fuel  expended  during  the  maneuver  from 
the  initial  circular  orbit  to  the  final  circular  orbit.  A  formulation  of  the  problem  is 
established  which  is  conducive  to  determining  the  overall  optimal  solution. 

A.  THE  PROBLEM 

The  overall  optimization  problem  for  the  orbital  plane  change  is  an  optimization 
problem  attempting  to  achieve  the  maximum  orbital  plane  change  using  the  available 
fuel.  The  total  fuel  to  carry  out  the  maneuver  consists  of  the  sum  of  the  propellant  used 
for  the  deorbit  impulse,  the  atmospheric  pass  and  the  reorbit  and  circularization  impulses. 
This  allows  many  possibilities  of  how  to  analyze  the  problem.  Since  the  exoatmospheric 
impulses  are  discrete  in  nature  they  could  be  thought  of  as  discrete  parameters.  They  can 
also  be  used  to  contribute  to  the  total  inclination  change  if  not  fired  completely  in  the 
plane  of  the  current  orbit.  It  has  been  shown  by  Miele  and  Lee  [Ref.  5]  that  by  applying 
the  exoatmospheric  impulses  slightly  out  of  plane,  improvements  in  performance  can  be 
achieved.  The  out  of  plane  yaw  angles  used  for  each  impulse  would  then  become 
additional  parameters  of  the  optimization  problem.  This  overall  problem  could  be 
treated  as  a  parameter  optimization  problem  attempting  to  maximize  the  total  inclination 
change  during  four  discrete  events,  the  deorbit  impulse,  the  atmospheric  pass  and  the 
reorbit  and  circularization  impulses.  In  this  case,  it  would  be  necessary  to  create  an 
algorithm  to  solve  for  the  optimal  control  history  of  an  atmospheric  pass  starting  with  the 
initial  atmospheric  entry  conditions  dictated  by  the  deorbit  impulse  and  terminating  at  the 
atmospheric  exit  conditions  dictated  by  the  reorbit  and  circularization  impulses.  The 
theory  which  will  be  considered  for  the  optimization  of  the  atmospheric  pass  will  be 
discussed  in  later  chapters.  The  remainder  of  this  chapter  will  be  devoted  to  the 
modeling  and  formulation  issues  present  in  the  overall  problem. 

(Note:  All  variable  definitions  are  stated  in  the  List  of  Symbols  at  the  front  of 

this  thesis) 
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B.  EXOATMOSPHERIC  PORTION 

The  three  impulses  which  are  performed  outside  the  atmosphere  would  determine 
the  entry  and  exit  conditions  for  the  atmospheric  optimal  control  problem.  For  example, 
once  the  <7^^  and  AF^orW/are  chosen,  the  equations  for  conserving  energy  and  angular 
momentum  will  determine  the  entry  velocity  and  flight  path  angle.  In  like  manner,  once 
the  AVreorhil  ,AFrirc  and  their  associated  <r‘s  are  chosen,  the  atmospheric  Vail  and 
will  be  fixed.  The  inclination  change  from  AFreorto  can  be  found  looking  at  the  vector 
velocities  as: 


=  tan' 


-if  ^Vdeorh„ 
f  y orbit  -&V<korb«  COS(^orto) 


The  other  inclination  changes  associated  with  the  impulses  can  be  found  in  like 


manner. 


Other  known  entry  and  exit  boundary  conditions  are  the  entry  and  exit  radii  which 
are  both  equal  to  .  The  entry  0,  <j>  and  \]/  angles  are  also  fixed  from  the  Keplerian 
elliptical  deorbit  at  the  radius  .  Also  the  entry  and  exit  masses  will  be  fixed  by  the 
fuel  expended  during  the  impulses.  Thus,  the  entire  entry  state  vector  is  fixed  by 
choosing  R^,  and  AV^,,.  Also,  a  target  state  space  is  set  since  m^,,  F^,, 

7 exi, and  Rexu  are  specified  by  the  exoatmospheric  parameters.  This  leaves  three  exit 
states  free  for  the  algorithm  to  use  in  the  optimization  process. 

C.  COORDINATE  SYSTEM  AND  EQUATIONS  OF  MOTION 

The  coordinate  system  and  equations  of  motion  used  here  are  for  a  lifting  body  in 
an  inverse  square  gravity  field.  These  equations  are  derived  by  Vinh  [Ref.  14],  and 
assume  a  spherical  central  body  and  a  non-rotating  atmosphere.  The  vehicle  is  a  point 
mass  experiencing  aerodynamic  forces. 

The  state  vector  is: 


[r  v  y  0  <|)  V)/  m]T 

The  control  vector  is: 


[a  5  a  T]  T 
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The  thrust  vector  is  fixed  along  the  reference  body  axis  from  which  a  is 
measured.  The  yaw  angle  is  only  used  as  a  control  in  the  exoatmospheric  portion  of  the 
maneuver  and  is  assumed  zero  during  atmospheric  flight. 

The  gravitational  field  is  : 

g  r2  (2) 


This  system  is  depicted  in  Figure  1  below. 


dv__T  cos(  a.)  -  D(  a.)  jul  sin(y) 

dt  m  r2  (4) 
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(5) 


vdj  (L(a)-f  7  sin(a))  cos(8)  ficos(y)  v2cos(y) 

dt  m  2  +  r 

r 


dQ  v  cos(y)  cos(q;) 
dt  rcos(<|)) 


d§_  v  cos(y)  sin(i|/) 
dt  r 


-  —  _  (k(a) +  Tsin(g.))  sin(S)  v  cos(y)  cos(ip)tan(<|)) 

dt  m  cos(y)  r 


L(a)  and  D(a)  are  computed  as: 


D(a)=-pv2AQ,(a) 


L(a)=-pv2AQ(a) 


From  the  rocket  equation  we  have  the  mass  state  derivative: 

dm  T 
dt  IspgQ 

The  orbital  inclination  at  any  time  is  given  by: 

cos(i)=cos(<j))cos(v|/) 


The  heat  rate  is  determined  by  using  Chapman’s  equation: 

^=Cp°V15 
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The  heating  rate  used  in  this  thesis  will  be:  100  -rr~ —  or  1. 13  xlO6  This 

ft  -  sec  m 

specific  value  was  used  here  for  comparisons  to  be  made  later. 

Density  is  modeled  as  locally  exponential  for  analytical  purposes,  that  is: 

p=P0<rAr~ro)  (14) 


D.  THE  PLANET  MODEL 

The  planet  model  used  throughout  this  section  is  a  spherical  earth  with  a  non¬ 
rotating  locally  exponential  atmosphere.  The  1962  U.S.  standard  atmosphere  is  used  in 
the  numerical  analysis  with  the  local  p  derived  from  that  data. 

E.  THE  VEHICLE  MODEL 

The  vehicle  used  throughout  this  thesis  is  the  MRRV  which  has  the  following 
characteristics.  [Refs.  2,4] 


3.96  m 


Figure  2:  MRRV  Diagram  from  [Ref.  2] 


Mass  (fully  fueled):  4899  kg 
Fuel  capacity:  2615  kg 
Aerodynamic  reference  area:  11.7  m2 

Engine  characteristics  and  Isp:  Three  Marquardt  R-40-B  rocket  motors  make  up 
the  propulsion  system,  providing  a  total  of  14679  N  of  thrust  at  a  specific  impulse  of  300 
sec. 
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Aerodynamic  Coefficients: 

The  aerodynamic  data  for  the  MRRV  comes  from  an  analytic  computer 
extrapolation  of  wind  tunnel  data  taken  by  the  Air  Force  Flight  Dynamics  Laboratory 
(AFFDL).  The  actual  wind  tunnel  tests  were  performed  on  a  scale  model  of  the  MRRV 
at  transonic  Mach  numbers  as  well  as  Mach  3  and  Mach  5  [Ref.  8],  This  original  data  is 
available  in  tabular  form  either  from  AFFDL  or  the  Naval  Postgraduate  School  library. 
1’his  data  was  then  adjusted  to  take  into  account  the  increase  in  Mach  number  and 
altitude  e  rerienced  in  the  orbital  maneuver  environment.  The  analytic  projections  are 
tabulated  in  terms  of  polynomial  fits  for  axial  and  normal  force  coefficients  in  [Ref.  4]. 
Although  this  data  is  given  for  various  Mach  numbers  the  only  data  used  in  this  thesis  is 
that  for  Mach  20+  since  the  entire  maneuver  takes  place  in  that  range.  A  typographical 
error  was  noted  in  [Ref.  4]  in  the  exponent  of  one  of  the  polynomial  coefficients.  The 
corrected  fit  data  is  listed  below.  The  axial  force  coefficient  is  composed  of  a  skin 
friction  term  and  a  pressure  term  that  is, 

CA=CAsF+CAfR  (15) 

where 

CAsF=A1h2  +  B1h  +  C1  (16) 

h=  210000  ft  was  used  as  a  representative  altitude,  giving  a  value  of  .0105  for 
the  skin  friction  term. 

CApR=A2a2+B2a  +  C2  (17) 


CN  —  A3a2  +  B3a  +  C3 
a  in  degrees,  0°  <  a  <  40° . 


(18) 
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The  values  used  here  are  given  below: 


subscript 

A 

B 

C 

1 

1.38E-12 

LO 

oo 

.0297 

2  (0°  <  a  <  16°) 

1.33E-5 

-8.19E-4 

.0243 

2(16°  < a <40°) 

-6.25E-6 

3.58E-4 

.0105 

3 

4.38E-4 

6.50E-3 

-.035 

The  one  piece  of  data  which  was  noted  as  an  error  in  [Ref.  4]  was  the  exponent 

for  Bs. 

The  data  from  these  polynomials  was  then  transformed  into  lift  and  drag 
coefficients  using  the  following  equations: 

cd  =  Ca  cos(a)  +  CN  sin(a)  (19) 

CL  =  CN  cos(a)  -  CA  sin(a)  (20) 

Once  transformed  a  cubic  polynomial  was  then  fit  to  the  lift  and  drag  coefficients, 
using  a  least  squares  fit,  yielding  the  following  polynomials: 

CD  =  (4.9852E-6)a3  +(2.2429E-4)a2  -(2.4732E-3)a  +  (3.6994E-2)  (21) 

CL  =  (-5.8340E-6)a3  +(5.6214E-4)a2  +  (5.0475E-3)a-(3.3418E-2)  (22) 

For  a  in  degrees,  0°  <  a  <  40° . 

These  are  depicted  below  as  well  as  the  lift  to  drag  ratio  for  use  as  the 
aerodynamic  model  for  this  thesis.  These  plots  are  of  the  curve  fit  data,  but  are 
indistinguishable  from  the  plotted  data  on  the  same  scale.  The  errors  in  the  fit  data  from 
the  data  in  [Ref.  4]  are  shown  in  Figure  6  through  Figure  8.  Figure  6  shows  the  percent 
error  in  the  lift  coefficient  fit  data  for  a  >  10  degrees.  The  percent  errors  are  larger 
below  10  degrees  because  the  lift  coefficient  goes  through  zero.  Therefore,  the  absolute 
error  is  shown  for  a  <  10  degrees  in  Figure  7.  The  percent  error  in  the  drag  coefficients 
data  is  shown  in  Figure  8. 


ll 


-3 


Figure  7:  Absolute  Error  in  Lift  Coefficient  vs.  Angle  of  Attack 


Figure  8:  Error  in  Drag  Coefficient  (%)  vs.  Angle  of  Attack 

In  order  to  solve  the  above  problem,  many  approaches  may  be  taken.  There  are 
various  direct  methods  and  some  “indirect  methods”  of  performing  this  maximization 
of  the  final  inclination.  The  next  chapter  will  consider  the  theory  of  an  indirect  approach 
and  the  application  of  the  indirect  method  to  solve  the  above  problem.  Then  some  direct 
method  analysis  will  be  considered  to  explore  what  is  learned  through  the  indirect  theory. 


14 


III.  THEORETICAL  ANALYSIS  (INDIRECT  METHOD) 


The  Restricted  Maximum  Principle  refers  to  that  version  of  the  Pontryagin’s 
Maximum  Principle  (PMP)  which  incorporates  state  variable  constraints.  This  material 
is  dealt  with  in  many  optimal  control  texts.  However,  the  notation  and  the  discussion  are 
such  that  it  is  necessary  to  review  multiple  texts  to  obtain  a  good  understanding  of  how  to 
formulate  this  problem.  This  chapter  will  establish  the  notation  and  definitions  to  be 
used  in  this  thesis  with  regard  to  this  Restricted  Maximum  Principle  (RMP)  which  is  a 
combination  from  various  sources.  [Refs.  9-11]  are  the  principal  texts  for  this 
formulation.  Also,  this  discussion  is  not  intended  to  serve  as  a  general  theoretical 
presentation  of  the  RMP  but  rather  as  background  for  applying  the  principle  to  this 
constrained  plane  change  optimization.  The  theorems  are  presented  without  proof.  See 
the  above  texts  for  more  detailed  discussions. 

A.  GENERAL  PROBLEM  FORMULATION 
Consider  the  autonomous  control  problem  given: 

a)  x  =  f(x,  u),  x  e  R",  ueR1  with  f  continuously  differentiable. 

b)  the  initial  point  x0  and  the  final  target  in  the  form  of  a  vector  valued  function 
P(x(f/)>//)eRm such  that  F=0  implies  that  the  target  is  hit, 

c)  the  class  of  admissible  controls,  which  is  taken  to  be  all  piecewise  continuous 
functions  u,  with  u(/)  eU.UcR5 

d)  the  state  constraint  which  requires  x(r)  to  lie  in  a  closed  region  B  of  the  state 
space  R”  of  the  form 

B  =  {  x:G(x)<  0}  (23) 

for  some  given  scalar-valued  function  G,  having  continuous  second  partial 
derivatives,  and  for  which 


VXG 


r  dQ 


dGr 


does  not  vanish  on  the  boundary  of  B,  that  is,  on  {  x:  G(x)  =  0} . 


(24) 
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Find:  The  control  histoiy  iT(/),  0  <  t  <  tf ,  that  minimizes  the  cost  functional 

'/ 

J(u>=  ),tf\ +  J/0(x(r),  u {t))dt  (25) 

0 

(Note  that  this  statement  of  the  theory  uses  a  Bolza  performance  index  for  generality. 

For  the  orbital  plane  change  problem,  /0  is  zero.) 

B.  THE  BASIC  MAXIMUM  PRINCIPLE 

The  Basic  Maximum  Principle  (BMP)  covers  the  above  problem  in  the  absence  of 
the  state  constraint  d).  This  principle  will  be  stated  first  and  then  modified  to  form  the 
Restricted  Maximum  Principle  (RMP). 

Define  for  X  e  R” 

H(x,  u,  X  )  =  -  /0(x,  u)  +  X  •  f (x,  u)  (26) 

(H  is  called  the  Hamiltonian)  and 

M(x,  X  )=  maxuelJ  H(  x,  u,  X  )  (27) 

Then  the  principle  can  be  stated  as:  [Ref.  Pontryagin  see  Knowles  page  57  [4]] 

Suppose  u  is  an  optimal  control  for  the  above  problem  and  x*  is  the 
corresponding  trajectory.  Then  there  exists  a  non-vanishing  function  X’(t) ,  hereafter 
referred  to  as  costates,  such  that 


(i)  x*  =  f(x\u*) 

(28) 

...  •.  d  H 

00  X  =  - 

0  X 

(29) 

(iii)  H(x\u*,X')  =  M(x\X*)=  maxueU  H(x\u,X*) 

(30) 
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The  boundary  conditions  to  be  satisfied  by  the  above  In  dimensional  system  of 
first  order  ODE  s  are  given  by  x0,  F=0,  and  the  terminal  transversality  conditions: 


d<S>  (d  F 


d  x  +  x 


d<&  (dt 


'='/ 

C.  THE  RESTRICTED  MAXIMUM  PRINCIPLE 

Now,  restrict  the  control  problem  further  applying  d)  from  above.  In  particular, 
suppose  the  set  of  admissible  controls  U  is  the  set  of  all  u  e  Rsfor  which: 

b,(u)<0,  b2(u) <  0, ...,  bt.(u) < 0,  (33) 


for  given  continuously  differentiable  functions  bj ,  b2>  ...  ,bA.in  particular 
U=[-l,l]  could  be  written 

{m  eR:«2  -1  <  0},  (34) 

that  is 

b,(«)=M2-l,  as R.  (35) 

If  we  define 

p(x,u)=VxG(x)-f(x,u),  (36) 

then  once  a  trajectory  x(t)  with  control  u(r),  hits  the  boundary  of  B,  G(x>=0,  a 
necessary  and  sufficient  condition  for  it  to  remain  there  for  all  later  t  is  that 

p(x(r),u(0)  =  0.  (37) 

This  asserts  just  that  the  velocity  of  a  point  moving  along  the  trajectory  is  tangent 
to  the  boundary  at  time  t .  Two  further  assumptions  that  will  be  required  later  are  that 

Vup(x,u)*0,  (38) 

and  that 

VuP(x,u),  Vub,(u),...,  Vnbfc(u)  (39) 

be  linearly  independent,  on  the  boundary  of  B,  along  an  optimal  trajectory  x*(/),u*(r). 
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This  ensures  that 


k  +  l<s  (40) 

The  strategy  for  solving  state  constrained  problems  is  to  break  them  into  two  sub¬ 
problems.  Whenever  the  optimum  trajectory  x*(/)lies  in  the  interior  of  B,  G(x’(/))<0, 
the  state  constraints  are  superfluous  (non-binding),  and  the  BMP  is  applied.  However,  as 
soon  as  x  (t)  hits  the  boundary,  G(  x*(r)  )=0,  the  RMP  is  applied.  In  this  way,  as  long  as 

x*  consists  of  a  finite  number  of  pieces  of  these  two  types,  there  is  hope  to  construct  it. 
The  essential  point  of  the  RMP  is  that  once  the  boundary  is  contacted,  the  control  power 
must  be  limited  to  prevent  the  optimum  trajectory  from  leaving  B.  From  the  above 
remarks,  this  requires  that  the  optimal  control  u*(r)  not  only  take  on  values  in  U,  but 
also  that 

p(x*(0,«*(0)  =  C)  (41) 

Define  also 

L=  H-JX-b.-rip,  (' 4 2) 

(L  is  called  the  Lagrangian  where  ri  (r)sK,  (/),.,  kk(0,  are  Lagrange  multipliers) 

The  RMP  takes  into  account  the  added  control  restraint  of  Equation  (41)  and  can 
be  stated  as: 

Let  x*(r),  t0  <f  <r„  be  an  optimal  trajectory  with  control  u*(r),  and  suppose 
that  x’(/)  lies  entirely  on  the  boundary  of  B,  t0  <  t  <  tx ,  and  Equations  (38)-(41)  are 
satisfied.  Then  there  exists  a  continuous  vector  function  X(t)  and  a  piecewise  smooth 
scalar  function  r|(r),  r0  <  t  <  tx  such  that 

(i)  x‘  =  f(x\u*),  (43) 

(ii)  1*  =  -VxH(x\u\*,)+r|(0  Vxp(x\u*)  ^ 
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(iii)  H(x\u\A.)  =  M(x',X)=  maxuelI  H(x\u,X)  (45) 

for  t0<t<  r„  where  the  maximum  is  taken  subject  to  uell,  that  is.  Equations 

(33)  and  (37).  Since  this  is  really  just  a  constrained  maximization.  Equation  (45) 
can  be  alternatively  written 

(iv)  V„L-0=V.H(x',u',X‘)-ri(f)  V„p(iV)-|;K,b,(H-)  (46) 

t0<t  <tx,  at  least  at  points  of  differentiability  of  q. 

(v)  A,(r0)  ^  0  (47| 

and  X(t0 )  is  tangent  to  the  boundary  of  B  at  x(t). 

(vi)  Whenever  q  is  differentiable  q(/)Vg(x*(0)  is  directed  toward  the  interior  of 
B  or  is  zero. 

From  the  assumptions  of  Equations  (38)-(41),  (iv)  can  be  discussed  in  more  detail.  At 
any  point  of  continuity  of  q,  for  given  x*(r),u*(r),  and  k(t),  (iv)  is  a  system  of  s 
equations  in  k  +  \  unknowns  q  (O.k/O.-.k^O,.  Assumptions  (38)  and  (39) 
guarantee  that  the  matrix  of  the  system  has  rank  k  + 1 ,  and  so  the  unknowns  can  be 
solved  for;  in  particular  q  can  be  solved  for  as 

q  ( t )  =  X(t)  •  a  (t),  r0  <  r  <  r, ,  (48) 

where  a  is  a  piecewise  smooth  function.  An  important  point  here  is  that  with  a  slate 
constraint  in  the  problem  k  + 1  must  be  less  than  or  equal  to  5  in  order  to  solve  the 
system,  limiting  the  number  of  control  constraint  functions  allowed. 

The  boundary  conditions  for  the  differential  system  in  the  RMP  not  changed  by 
the  state  constraint.  What  is  changed  here  is  the  dynamics  of  the  system  that  must  meet 
the  boundary  conditions. 

D.  JUMP  CONDITIONS 

As  pointed  out  above,  the  strategy  for  constructing  the  optimal  trajectory  is  to  find 
a  finite  sequence  of  alternating  subarcs  which  are  either  not  affected  by  the  state 
constraint  or  lie  along  the  boundary  of  the  constraint.  This  will  require  the  subarcs  to  be 
connected.  At  these  connecting  points  additional  conditions  are  imposed  on  the  optimal 
trajectory.  These  additional  conditions  are  called  ‘jump  conditions’  because  they 
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potentially  force  the  costates  to  have  a  jump  discontinuity  upon  exit  from  the  boundary  of 
B.  The  costates  are  assumed  to  be  continuous  upon  entering  the  constrained  arc  and 
then  a  modification  must  be  made  to  the  costate  values  on  exiting  the  constrained  arc.  If 
{a  is  the  time  the  constraint  is  encountered  and  /Ais  the  time  that  the  constrained  arc  is 
terminated  the  jump  conditions  are  as  follows: 


M>.-  o)  =  a.(/„+0). 

(49) 

Hh  +  0)  =  Mtt  -  0)- iiMVMMt,)) 

(50) 

E.  THE  COMPLETE  FORMULATION 

In  order  to  generate  a  complete  trajectoiy  from  t  =  0  to  /  =  tf ,  the  equations 

would  be  propagated  as  an  unconstrained  arc  applying  the  Basic  Maximum  Principle 
until  the  boundary  was  reached.  At  that  tune  the  equations  would  be  propagated  using 
the  RMP  and  upon  leaving  the  boundary  of  B  the  jump  conditions  would  be  applied  to 
the  costates  and  the  equations  would  continue  applying  the  basic  principle  again,  etc.  By 
this  method  a  series  of  alternating  arcs  is  generated  which  when  satisfying  all  of  the 
boundary  conditions  produces  the  optimal  trajectoiy. 

A  note  about  the  Lagrange  multipliers  mentioned  in  the  RMP.  The  Lagrange 
multipliers  have  the  property  that  they  are  zero  when  the  associated  constraint  is  not 
binding  and  non-zero  when  the  constraint  function  is  equal  to  zero.  In  the  case  of  q,  q=0 
implies  that  G(x'(7))<0  which  is  when  the  basic  principle  is  applied  and  q  is  non-zero 
when  G(x  (/ )  )=0  and  the  RMP  is  applied.  The  k’s  have  this  property  relative  to  their 

associated  q  as  well  as  the  restriction  that  they  be  positive  when  their  associated  q=0. 
This  can  be  written  as 


-=5 

XJ 

li 

o 

(51) 

K.b,  =  0 

(52) 

K;  >  0 

(53) 

b,.  <0 

(54) 

This  property  is  sometimes  referred  to  as  complimentary  slackness. 
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So,  at  this  point  there  is  a  2  n  dimensional  system  of  first  order  ODE’s.  The 
2/7  +  1  boundary  conditions  are  provided  by  n  initial  conditions  x0,  m  target  conditions 

from  F(x(tf),tf)  e  R“ ,  and  the  remaining  n  +  \-m  conditions  come  from  the  terminal 

transversality  conditions.  The  costate  equations  have  jump  discontinuities  as  they  exit 
the  constrained  arcs  which  are  specified  by  the  jump  conditions.  The  difficult  task  at  this 
point,  as  is  almost  always  the  case  when  using  indirect  methods,  is  solving  this  two  point 
boundary  value  problem  (TPBVP). 

F.  APPLICATION  OF  THE  RESTRICTED  MAXIMUM  PRINCIPLE 

Now  the  above  theory  will  be  applied  to  the  atmospheric  pass  of  the  problem  at 
hand.  The  state  vector  is  modified  to 

x  =  [/•  v  y  <p  ij j  m]T  (55) 

i.e.,  the  0  state  is  left  out.  This  is  because  it  turns  out  that  the  0-costate  Xa  is 
identically  zero  for  the  entire  trajectory  during  both  the  constrained  and  unconstrained 
portions.  This  comes  from  the  fact  that  the  0-costate  derivative  is  identically  zero  and 
the  terminal  transversality  condition  gives  the  0-costate  a  value  of  zero  at  the  terminal 
time  making  this  costate  identically  zero.  Since  the  0— costate  is  zero  and  0  does  not 
appear  in  any  other  state  derivatives,  it  need  not  be  calculated  in  determining  the  optimal 
control  for  this  problem.  It  can  be  computed  after  the  final  iteration  from  the  other  states 
if  needed  to  visualize  the  trajectory.  This  reduces  the  problem  from  14  equations  to  12 
and  significantly  simplifies  the  remaining  costate  derivatives.  For  applying  the  RMP  we 
have  the  following: 

The  cost  functional: 

J  =cos(i(tf  ))=^>[x(tf),tf]  =cos(<K  tf  »cos(vf i(tf ))  (56) 

Note  that  f0  -  0  in  this  case. 

The  costate  vector: 

X  =  IKK  \  \  K  Kf  (57) 

The  control  vector: 

u  =  [T,ct,8f  (58) 
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The  admissible  controls: 

U={u:0<T<T_,0<a<amax}  (59) 

The  control  constraint  function: 

The  constraints  on  T  and  a  could  be  represented  by  two  functions: 

brf-Tl_  (60) 

and 

b2=a2-aamax  (61) 

However,  in  this  case,  as  will  be  seen  later  in  Equation  (92),  the  derivative  of  L 

with  respect  to  8  does  not  contain  any  of  the  Lagrange  multipliers.  This  effectively 
reduces  the  number  of  equations  in  the  system  to  be  solved  to  2.  This  means  that  we  are 
only  allowed  two  Lagrange  multipliers  if  the  system  is  to  be  solved.  Thus,  we  must  find 
one  function  which  encompasses  all  of  the  constraints  on  T  and  a.  Suppose  we  define  b 
as 


-1 


(62) 


where  z  is  large,  positive,  and  even.  This  exponent  causes  the  function  to  take 
on  a  value  of  almost  exactly  -1  whenever  0  <  T  <  T  max  and  0  <  a  <  a^,  and  a  value  of 
almost  exactly  0  when  either  control  is  at  its  boundary  with  the  other  in  the  interior.  The 
value  of  b  is  very  large  and  positive  when  any  control  is  outside  the  allowed  region.  By 
choosing  z  arbitrarily  large  (e.g.  z  =  200)  this  function  can  be  made  to  come  arbitrarily 
close  to  the  comers  of  the  region.  The  following  figures  will  demonstrate  the  effect  of 
increasing  z .  For  each  value  of  z  two  plots  are  shown.  The  first  is  a  3-d  view  of  the 
function  shape.  The  second  is  a  2-d  view  for  the  angle  of  attack  side  of  the  function. 
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control  constraint  function  b 


The  state  derivatives: 


x  =  f(x,u) 

As  described  in  the  equations  of  motion  in  the  previous  chapter. 
The  state  constraint: 

?=cPv5<c 

The  state  constraint  function  ( G  <  0  =>  q  <  gmax ): 

where  N=.5,  M=3. 15  and  qmax  is  a  positive  constant. 
The  boundary  function  (p  =  Q  =>  q  ~  const ): 

/>  =  VxG-f(x,u) 


p  =  rcos(a)-D(a)-msin(y)  -^-  +  £v2 

l  r2 


where 


s  3.15 


The  Lagrangian:  (  L=  X- f  (x,  u)  -  Kb  -  rjp): 


l  =  X  vsm(r)  +  A.  fl£25<£)zSte)_M5(I) 

r  V  m  2 


^  (L(a)  +  rsin(q))  cos(S)  p.cos(y)  y_  cos(y) 


X  v  cos(y)  sin(qf) 


^  (L(q)  +  T  sin(q))  sin(S)  vz  cos(y)  cos(  qt )  tan(  $ ) 

vl _ m  cos(y)  r 


IspgO 


(  T  f  CL  ^  (  f 

2  J  ~ 1  +  2“ - -1  -I  Jcos(a)-'D(a)-msin(y)M~  +  ^v2 

A  max  )  l  max  J  J  V  l  r2 
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The  costate  derivatives :  ( X  - - ) 

d  x 


.  M-cos(y)  v  cos(y) 

dX  X  p.sin(y)  y  3  2  X  v  cos(y)  sin(u/) 

r  _  ~  y  '  V  r  ry<b 


X  vcos(y)cos(ij;)tan(<t)) 

„  T]  wsin(y)  jll 

2  3 
r  r 


(L(a)-f  7sm(g))cos(S)  jicos(y)  vzcos(y) 

dK  V°SW  r  »  3  +  r 

_=-if!h(y)-2J_+_ i - 5 - - - 

V 

V  cos(y)sin(q;)  X  cos(y)  cos(ij/)tan(<|>) 

- +  2  — - 

r  t 


^  (L(a)  +  7 sin(a.))  sin(S)  vz  cos(y)  cos(\|/)tan((|>) 

ty  y  m  cos(y)  r 


2  T\m  sin(y)  \  v 


f  .  ■>  \ 

jisin(y)  v  sin(y) 

^  M-cos(7)  y  2  r  X  vsin(y)sin(*|/) 

-ji—X  »<«(?)+-!— - i— ^ - i+jf - 

r  ^.2  v  r 

/  ~  N 

^  (L(«.)  +  rsin(a))  sin(S)  sin(y)  ^  sin(y)  cos(t|;)tan((|)) 

^  t _ m  cos(y)2  r  , 


,  .  [  LL  7 

-  T|  m  cos(y)  —  +  §  v 
;  r2 


dX^  A^vcos(y)cos(t|/)(l+tan(<j>)2) 
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^  v  cos(y)  cos(ip)  A,  v  cos(y)  sin(t|/)  tan(<|>) 


d%m  xv(rcos(a)-D(a))  *  (L(a)  +  rsin(a))cos(8)  %  (L(a)  +  rsin(a))sm(5) 


-T]sin(y)f-^-  +  £v2 


mr  cos(y)  v 


The  initial  state  vector  x0  is  completely  specified  as  discussed  in  the  overall 
problem  statement.  The  initial  states  are  all  determined  by  the  exoatmospheric  portion  of 
the  descent.  The  choice  of  cr(kortit  and  AV^oriil  will  specif/  the  atmospheric  entry  states. 
The  target  function: 

~r if/)' 

v  v  fJ  fJ  V  -v(t^) 

<76> 

mexit  ) 

These  conditions  are  specified  by  the  reorbit  and  circularization  impulses  and 
bank  angles  chosen. 

The  terminal  transversality  conditions:  ( X(tf  )=-  —  ^  +  f  )  v  I  ) 

\dx)  ) 

'='/ 

\(V)=V  i  (77) 

Xv(V)=V2  (?8) 

W)-V  3  (?9) 

Vy)=^*(y))ciXy))  (80) 

(8,) 


28 


The  above  equations  yield  two  non-linear  boundary  conditions  in  (80)  and  (81) 


The  terminal  Hamiltonian  value:  ( H {tf )  = 


) 


H(r/)  =  °  (83) 

It  should  be  noted  here  that  in  value  the  Hamiltonian  and  the  Lagrangian  are 
equal.  So,  this  is  equivalent  to  making  the  Lagrangian  zero  at  tf . 

Thus  there  are  13  boundary  conditions  for  the  12  equations  which  are  necessary 
to  be  satisfied  by  the  optimal  solution.  The  13*  is  necessary  because  the  terminal  time  is 
unknown.  It  is  possible  with  a  change  of  dependent  variable  and  the  addition  of  an 
additional  parameter,  namely  the  terminal  time,  to  cause  the  system  to  always  be 
integrated  between  0  and  1.  This  is  done  by  letting  x=t/tf ,  the  new  dependent  time 

variable  t,  is  then  the  variable  of  integration  and  all  12  of  the  differential  equations  have 
their  right  hand  sides  multiplied  by  tf .  This  in  essence  scales  all  the  equations  into  the 

0  to  1  time  frame.  This  is  a  good  idea  here  since  it  simplifies  the  boundary  value 
problem  by  fixing  the  terminal  time  at  1  at  the  cost  of  only  one  additional  parameter  in 
the  parameter  optimization  process  going  on  external  to  the  boundary  value  solver. 

The  jump  conditions  show  how  the  costate  variables  must  jump  when  exiting  the 
boundary  region  during  the  integration. 

The  jump  conditions:  (  X(tb  +  0)  =  X(tb  -  0)  -  p(/A  )VxG(x(tb ))  ) 


Xr(^  +  O)  =  Ar(^~O)+TlCp0pjVe^  ^  °^vM 

(84) 

\(f*+0)-\('4-0) - “ - 7 - 

(85) 

(86) 

(87) 
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w°)“\fr-°)  <88) 

x  (t.-ti)  (89) 

m\  b  )  m\  b  ) 

Here  again  it  can  be  seen  that  the  value  of  rj  is  needed  to  perform  the  integration. 
Also  note  that  the  only  costates  affected  are  those  whose  respective  states  appear  in  G. 
The  remaining  costates  are  continuos  at  both  entry  to  and  exit  from  a  constrained  arc. 

The  determination  of  the  value  of  r|  and  of  the  optimal  controls  go  hand  in  hand. 
The  necessary  conditions  on  the  controls  are: 


X  cos(a)  X  sin(a)  cos(8)  X  sin(a)  sin(S) 
v  y  v| / 

- 1 — i - _j — x - 

m  mv  m  cos(y)v 


—  1  z 


T  2  — 
max  T 


Tjcos(a) 


X  (L(a)  +  rsin(a))sin(8)  X  (L(a)  +  rsm(a))cos(S) 

_J _ +  _V _ 

mv  m  cos(y)  v 


_  X  -T  sin(a)  - 
d__  v\ _ 

da  m 


)JJ  ^(^Ua)j+rcoS(a)JcoS(8) 

+  m  v 


\  UaaI^aV  +  7 cos(a*J  sin(8) 

m  cos(y) v 


K  2 - 1  z 

a 

L  max  J 


a  2 - -1 

max  a 
V  max 


30 


From  Equation  (92)  a  relationship  for  the  optimal  bank  angle  can  be  derived: 


8=  arctan 


X  cos(y) 

k  y 


(94) 


This  relationship  holds  for  the  entire  trajectory,  both  on  constrained  and 
unconstrained  arcs. 


Equations  (91)  ,  (93),  (37)  and  the  complementary  slackness  property  given  in 
Equations  (51)  through  (54)  will  then  form  a  system  of  equations  in  the  unknowns 
needed  to  determine  the  optimal  values  of  T,  a,  k,  and  q.  A  difficulty  could  arise  if  k 
were  allowed  to  be  zero  for  any  finite  time.  The  difficulty  is  that  any  value  of  T  would 
satrsfy  these  conditions.  It  can  be  observed  that,  in  Equation  (91),  T  cannot  be  evaluated 
if  k=0.  When  the  second  derivative  of  L  with  respect  to  T  is  considered,  it  can  be  seen 
that  the  second  order  condition  for  maximizing  L  with  respect  to  T  can  only  be  satisfied 
when  k>0.  This  is  the  case  when  b=0. 


d  2L 
dT2 


-4k 


2  T 
T 

max 


-yv.) 


(95) 


Barring  the  possibility  of  k=0,  which  would  allow  any  value  of  T  to  be  optimal 
and  lead  to  a  so-called  singular  arc,  the  constrained  arcs  which  meet  these  conditions 
(i.e.,  b=0  and  p=0)  are  as  follows: 

(1)  T=0  with  a  satisfying  Equation  (67)  with  0  <  a  <amax. 

(2)  T=Tmax  with  a  satisfying  Equation  (67)  with  0  <  a  <  amax . 

(3)  a=0  with  T  satisfying  Equation  (67)  with  0  <  T  <  Tmax . 

(4)  ot=amax  with  T  satisfying  Equation  (67)  with  0<T<  Tmax . 

In  the  next  chapter  a  numerical  analysis  will  be  performed  to  show  how  this 
information  should  be  incorporated  into  direct  method  analysis  of  this  problem. 
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IV.  NUMERICAL  ANALYSIS  (DIRECT  METHOD) 


The  previous  chapters  have  given  some  idea  of  what  can  be  learned  by 
considering  the  indirect  methods  of  analysis  for  this  problem.  This  chapter  will  discuss 
a  direct  method  of  analysis  which  takes  into  account  the  conditions  which  must  be 
satisfied  by  the  optimal  trajectory.  This  narrows  the  possibilities  for  the  direct  method 
search.  Only  those  arcs  which  are  allowed  in  the  indirect  method  will  be  considered  here 
in  the  Aerobang  trajectory.  The  Aerocruise  analysis  is  done  for  comparison  to  results 
obtained  in  [Ref.  3]  which  did  not  take  into  account  the  existence  of  some  of  the 
previously  mentioned  arcs.  Many  of  the  parameters  of  the  problem  are  chosen  to  allow 
this  comparison.  This  chapter  will  also  discuss  how  this  direct  method  of  analysis  can 
aid  in  the  solution  of  the  indirect  problem.  Solving  the  indirect  problem  is  the  only  way 
to  actually  claim  optimality.  The  drawback  to  the  direct  analysis  is  that  no  matter  how 
good  the  trajectory  is,  there  is  no  way  to  claim  optimality  unless  that  trajectory  can  be 
shown  to  solve  the  necessary  conditions  discussed  previously. 

A.  POST 

The  program  which  was  used  for  the  direct  analysis  is  the  Program  to  Optimize 
Simulated  Trajectories  (POST).  POST  is  a  general  purpose  FORTRAN  program  for 
simulating  and  improving  point  mass  trajectories  of  aerospace  vehicles.  The  program 
can  be  used  to  examine  a  wide  variety  of  performance  and  mission  analysis  problems  for 
atmospheric  and  orbital  vehicles. 

The  basic  simulation  flexibility  is  achieved  by  decomposing  the  trajectory  into  a 
logical  sequence  of  simulation  segments.  These  trajectory  segments,  referred  to  as 
phases,  enable  the  trajectory  analyst  to  model  both  the  physical  and  non-physical  aspects 
of  the  simulation  accurately  and  efficiently.  By  segmenting  the  mission  into  phases,  each 
phase  can  be  modeled  and  simulated  in  a  manner  most  appropriate  to  that  particular 
flight  regime.  For  example,  the  planet  model,  the  vehicle  model,  and  simulation  options 
can  be  changed  in  any  phase  to  be  compatible  with  the  level  of  detail  required  in  that 
phase.  [Ref.  12] 
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B.  TRAJECTORY  CHARACTERIZATION 

The  POST  program  takes  a  single  formatted  input  file.  The  specifics  of  the 
formatting  requirements  are  discussed  in  [Ref.  13j.  It  is  sufficient  here  to  say  that,  as 
stated  above,  for  POST  to  perform  its  optimization,  the  desired  trajectory  must  be  broken 
down  into  phases  which  must  occur  in  order  and  as  such  constitute  an  a  priori  knowledge 
of  the  type  of  trajectory  POST  will  attempt  to  optimize.  Also  included  in  the  POST  input 
file  are,  the  control  parameters,  any  constraints  which  must  be  met  for  any  feasible 
trajectory,  and  the  values  of  the  control  parameters  to  be  used  as  the  initial  iterate.  The 
term  control  parameter  in  this  context  is  used  to  refer  to  those  parameters  which  POST  is 
allowed  to  perturb  in  its  search.  Two  different  models  were  examined  in  this  thesis  for 
the  characterization  of  the  orbital  plane  change  maneuver,  one  for  the  Aerobang 
maneuver  and  the  other  for  an  Aerocruise  maneuver.  The  phase  structure  will  be 
discussed  below  for  each  model  with  the  phases  listed  in  sequential  order.  The 
discussion  will  include  the  control  parameters  used  in  each  phase.  In  both  cases,  the 
target  was  to  have  zero  propellant  remaining  upon  trajectory  termination.  These  models 
represent  the  final  product  and,  as  will  be  discussed  later,  require  large  amounts  of  user 
interface  to  achieve.  Figure  15  shows  graphically  what  is  required  for  the  trajectory 
breakdown  with  the  numbers  representing  phases  of  the  trajectory. 


Earth 


Figure  15:  Example  Trajectory  Characterization  (Aerobang) 
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(1)  Aerobang  phases 

Phase  1  Initial  Orbit 

The  initial  orbit  for  both  models  is  a  circular  equatorial  orbit  at  185  km.  This 
initial  altitude  was  chosen  for  comparison  to  [Ref.  3],  In  this  phase  the  vehicle 
parameters  are  also  input. 

Phase  2  Deorbit  Impulse 

This  phase  begins  5  seconds  after  the  orbit  is  initialized.  An  impulse  in  the 
direction  of  the  negative  velocity  vector  is  applied  which  causes  the  orbit  to 
become  elliptical  with  apogee  at  185  km  and  perigee  less  than  111  km,  the 
altitude  at  which  the  atmosphere  will  be  ‘turned  on’.  Again,  these  altitudes  were 
chosen  for  comparison  to  [Ref.  3],  The  motion  is  determined  with  no  atmosphere 
during  this  phase.  The  magnitude  of  the  deorbit  impulse  is  a  control  parameter 
which  POST  has  at  its  disposal  for  its  optimization. 

Phase  3  Atmospheric  Entry  and  Glide  to  Heating  Rate  Boundary 
This  phase  begins  when  the  altitude  of  the  vehicle  reaches  1 1 1  km.  At  this  point 
the  equations  of  motion  begin  to  take  into  account  atmospheric  forces.  The  1962 
U  S.  standard  atmosphere  is  applied.  This  atmospheric  model  was  chosen  for 
consistency  with  [Ref.  3],  Thrust  is  zero  in  this  phase.  The  bank  angle  and  angle 
of  attack  are  used  as  control  parameters  by  POST.  These  two  controls  are 
represented  by  a  linear  polynomial  in  this  phase.  POST  implements  this  by  using 
the  constant  term  and  the  linear  term  of  these  two  polynomials  as  control 
parameters. 

Phase  4  Zero  Thrust  Constant  Heating  Rate  Arc 

This  phase  begins  when  the  heating  rate  reaches  a  value  equal  to  the  maximum 
allowed  heating  rate.  During  this  phase,  the  thrust  is  still  zero,  the  angle  of  attack 
is  modulated  according  to  Equation  (37)  to  keep  the  heating  rate  constant.  Note 
that  this  arc  can  only  be  accomplished  during  a  time  while  the  flight  path  angle  is 
negative  which  is  why  this  phase  is  assumed  to  be  the  beginning  of  the  heat  rate 
control  arc.  The  angle  of  bank  is  a  linear  polynomial  as  above  with  the  constant 
and  linear  terms  used  as  control  parameters.  In  order  for  POST  to  correctly 
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simulate  the  constant  heating  rate  axe  it  was  necessary  to  allow  POST  to  use 
angles  of  attack  greater  than  40  degrees  during  the  iterative  process  as  long  as 
angle  of  attack  was  less  than  40  degrees  for  the  entire  final  trajectory.  An 
additional  parameter  used  in  the  phase  is  the  duration  of  the  arc  itself.  Thus, 
POST  is  allowed  to  vary  the  length  of  this  phase. 

Phase  5  Max  Angle  of  Attack  Transition  at  Constant  Heating  Rate 
This  phase  was  included  as  a  result  of  the  transition  between  the  zero  thrust  and 
maximum  thrust  arcs.  The  thrust  is  throttled  during  this  phase.  As  discussed  in 
the  theoretical  treatment  of  the  problem,  this  arc  must  have  a=0  or  40  degrees. 
Since  a  is  less  than  90  degrees.  Equation  (37)  has  a  jump  discontinuity  in  a  if 
there  is  a  jump  in  T.  This  is  the  case  in  the  transition  from  zero  to  maximum 
thrust.  It  was  discovered  that  this  was  causing  a  to  exceed  40  degrees  at  this 
transition  in  the  final  trajectory.  This  phase  was  chosen  since  it  is  allowed  by  the 
indirect  analysis  discussed  previously  and  it  allows  the  vehicle  to  throttle  the 
thrust  and  build  up  gradually  to  maximum  thrust  while  not  exceeding  the 
maximum  angle  of  attack.  The  duration  of  this  phase  is  very  short  and  its 
addition  resulted  in  very  little  change  in  performance.  The  angle  of  attack  is  40 
degrees  and  constant  during  this  phase.  Bank  angle,  as  a  linear  polynomial,  is 
used  for  two  control  parameters  during  this  phase.  Thrust  is  modulated  to 
maintain  the  heating  rate  constant  until  maximum  thrust  is  reached. 

Phase  6  Max  Thrust  Constant  Heating  Rate  Arc 

This  phase  begins  when  the  thrust  increases  to  maximum  during  the  transition  in 
phase  5.  Thrust  is  maintained  at  maximum  and  a  is  modulated  to  maintain 
heating  rate  constant.  Bank  angle,  as  a  linear  polynomial,  is  used  for  two  control 
parameters  during  this  phase.  The  phase  duration  is  also  a  control  parameter  in 
this  phase. 

Phase  7  Glide  to  Atmosphere  Exit 

This  phase  begins  when  the  phase  duration  of  phase  6  is  reached.  The  goal  of  this 
phase  is  to  reach  the  top  of  the  atmosphere  at  111  km.  Thrust  is  zero  unless  it 
becomes  apparent  that  the  trajectory  will  not  make  it  to  the  top  of  the  atmosphere. 
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The  criteria  used  to  determine  this  is  the  flight  path  angle  equal  to  zero.  The 
phase  always  begins  with  a  positive  flight  path  angle  and  if  it  reaches  zero  the 
optional  phase  8  is  begun.  Angle  of  attack  and  bank  angle,  as  linear  polynomials, 
are  used  for  four  control  parameters  during  this  phase. 

Phase  8  (Optional)  Max  Thrust  Boost 

This  phase  begins  if  the  flight  path  angle  reaches  zero  during  phase  7  and  takes 
place  concurrently  with  phase  7.  The  rocket  engine  is  turned  on  and  the 
osculating  apogee  is  raised  to  185  km.  Angle  of  attack  and  bank  angle  are 
continued  as  controlled  in  phase  7.  In  POST,  this  event  is  a  sub-event  to  phase  7 
and  no  direct  control  is  introduced  here.  POST  continues  to  use  the  controls  it 
has  in  phase  7.  This  phase  terminates  when  the  apogee  altitude  reaches  185  km 
And  phase  7  continues  to  the  top  of  the  atmosphere. 

Phase  9  Atmospheric  Exit 

This  phase  takes  place  when  vehicle  altitude  reaches  1 1 1  km.  The  atmosphere  is 
turned  off  and  the  trajectory  continues  until  altitude  reaches  185  km.  No  control 
parameters  are  used  during  this  phase. 

Phase  10  Reorbit  Impulse 

This  phase  begins  when  the  flight  path  angle  reaches  zero  during  phase  9  and 
takes  place  concurrently  with  phase  9.  The  rocket  engine  is  turned  on  and  the 
osculating  apogee  is  raised  to  185  km.  This  is  necessary  because  the  drag  from 
the  atmosphere  causes  the  apogee  altitude  to  drop  even  after  the  boost  maneuver 
in  phase  8.  The  angle  of  attack  and  bank  angle  are  zeroed  and  the  impulse  is 
applied.  Since  POST  cannot  apply  an  impulse  whose  magnitude  depends  on  the 
current  orbital  parameters,  it  was  necessary  to  approximate  an  impulse.  This  is 
done  by  using  a  propulsive  rocket  firing  of  1000  times  the  magnitude  of  the 
normal  thrust.  This  allows  the  fuel  expended  during  this  impulse  to  be  taken  into 
account  and  the  firing  to  be  terminated  when  the  desired  orbital  parameters  are 
met.  This  firing  takes  a  very  short  time,  on  the  order  of  .1  seconds  and  closely 
models  a  pure  impulse.  In  any  case  this  firing  has  more  losses  than  an  impulse 
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and  is  conservative  in  terms  of  fuel  consumption.  No  control  parameters  are 
associated  with  this  phase. 

Phase  11  Circularization  Impulse 

This  phase  begins  when  altitude  reaches  185  km.  A  impulsive  thrust  is  applied  to 
circularize  the  orbit  at  185  km.  No  control  parameters  are  associated  with  this 
phase. 

Phase  12  Trajectory  Termination 

The  trajectory  is  terminated  as  soon  as  the  orbit  is  circularized.  POST  requires  a 
phase  to  exist  to  tell  it  that  the  trajectoiy  is  ended. 

(2)  Aerocruise  phases  ( Phases  1-3  are  identical  for  both  maneuvers.  Aerocruise 
phases  6-11  are  identical  to  Aerobang  phases  7-12.) 

Phase  1  Initial  Orbit 
See  Aerobang. 

Phase  2  Deorbit  Impulse 
See  Aerobang. 

Phase  3  Atmospheric  Entry  and  Glide  to  Flight  Leveling  Maneuver 
See  Aerobang. 

Phase  4  Flight  Leveling  Maneuver 

This  phase  begins  when  a  chosen  heating  rate  is  achieved.  This  is  one  of  the 
parameters  that  needed  to  be  adjusted  manually  to  target  the  heating  rate  desired 
for  the  cruise  arc.  Angle  of  attack  is  made  maximum  and  the  linear  term  of  the 
bank  angle  (bank  rate)  is  used  as  a  control  parameter  to  reduce  the  bank  angle  and 
cause  the  flight  path  angle  to  level  to  zero  for  a  steady  state  cruise. 

Phase  5  Steady  State  Cruise  at  Constant  Heating  Rate 

This  phase  begins  when  the  prescribed  flight  path  angle  is  reached.  Steady  cruise 
requires  flight  path  angle  to  be  identically  zero  for  the  entire  arc.  This  is  achieved 
by  initiating  the  arc  at  zero  flight  path  angle,  modulating  thrust  to  drag,  and 
modulating  bank  angle  to  maintain  y  equal  to  zero.  Since  POST  uses  input  as 
look-up  tables  in  which  it  uses  linear  interpolation  it  is  not  possible  to 
numerically  maintain  y  equal  to  zero.  This  results  in  small  changes  in  flight  path 
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angle  and  thus  altitude  and  density  in  the  simulated  runs.  A  slightly  positive 
flight  path  angle  was  used  to  initiate  this  phase  in  order  to  spread  the  numerical 
error  over  the  entire  arc.  The  variation  in  heating  rate  is  small,  about  1.5  %,  and 
average  heating  rate  on  this  arc  was  looked  at  which  is  conservative  in  terms  of 
trajectory  performance.  Angle  of  attack  is  assumed  constant  for  steady  state 
cruise  and  only  this  constant  term  and  the  phase  duration  are  used  as  control 
parameters  in  this  phase. 

Phase  6  Glide  to  Atmosphere  Exit 
See  Aerobang. 

Phase  7  Max  Thrust  Boost 
See  Aerobang. 

Phase  8  Atmospheric  Exit 
See  Aerobang. 

Phase  9  Reorbit  Impulse 
See  Aerobang. 

Phase  10  Circularization  Impulse 
See  Aerobang. 

Phase  11  Trajectory  Termination 
See  Aerobang. 

C.  THE  MAN  IN  THE  LOOP 

As  alluded  to  above,  working  with  POST  is  not  a  hands-off  experience.  The  user 
must  have  a  thorough  knowledge  of  the  problem  being  investigated  in  order  to  interpret 
the  results.  The  following  is  a  small  list  of  some  examples  of  the  type  of  intervention 
required  by  the  user. 

(1)  As  discussed  above,  the  problem  must  be  set  up  in  phases.  This  can  be  done 
in  multiple  ways.  In  [Ref.  1],  Nicholson’s  characterization  of  the  two 
maneuvers  is  quite  different  from  what  was  done  here.  Many  changes  were 
made  to  ‘clean  up’  the  final  trajectory  produced.  The  deorbit  and 
circularization  impulses  were  modeled  in  [Ref.  1]  as  propulsive  retro  firings  at 
180  degrees  angle  of  attack.  These  were  changed  to  pure  impulses.  For  the 
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Aerobang  problem,  the  angle  of  attack  required  to  maintain  the  heating  rate 
constant  is  not  always  less  than  maximum  angle  of  attack.  This  can  be  dealt 
with  in  multiple  ways.  [Ref.  1]  began  the  heat  rate  control  at  a  time  in  the 
trajectory  determined  by  previous  runs  and  limited  the  angle  of  attack  to 
maximum.  This  resulted  in  not  exactly  hitting  the  heating  rate  and  limited  the 
possible  trajectories  tested  for  the  optimization.  This  was  changed  to 
initiating  the  heat  control  phase  at  the  max  heating  rate  and  allowing  the  use 
of  higher  values  of  ct  than  40  degrees.  An  inequality  constraint  was  then 
placed  on  the  maximum  value  of  a.  This  allowed  a  larger  space  of  potential 
trajectories.  [Ref.  1]  did  not  include  the  zero  thrust  heat  rate  control  arc.  This 
was  added  to  again  increase  the  space  of  trajectories  available  for  the 
optimization.  It  was  assumed  in  [Ref.  1]  that  the  atmospheric  J3  was  a 
constant,  which  is  not  the  case,  resulting  in  non-constant  heating  rate  arcs. 
This  was  changed  to  use  the  data  from  the  1962  atmosphere  to  produce  a  table 
for  (3.  There  is  not  necessarily  a  “right  way”  or  a  “wrong  way”  to  set  up  ihe 
trajectory  but  the  results  can  easily  be  limited  by  the  users  characterization. 

(2)  The  results  of  POST’s  optimization  is  heavily  dependent  on  the  initial 
trajectory  as  well  as  the  overall  characterization.  For  instance,  if  the  angle  of 
attack  is  initially  decreasing  during  the  initial  glide  in  the  atmosphere,  POST 
searches  decreasing  angle  of  attack  trajectories,  while  the  opposite  is  true  if 
the  angle  of  attack  is  initially  chosen  to  be  increasing.  This  is  because  POST 
performs  a  local  maximization  of  the  optimization  variable.  As  a  result, 
POST’s  final  answer  may  very  well  be  a  local  extrema  and  the  POST  output 
may  indicate  that  the  problem  is  solved,  but  the  answer  may  clearly  not  be  the 
best.  For  this  reason,  it  was  necessary  to  try  various  starting  point  in  terms  of 
initial  guesses  to  allow  more  possibilities  to  be  searched.  Another  aspect  of 
the  initial  guess  is  that  sometimes  POST’s  search  may  lead  it  to  a  region  were 
it  cannot  target  the  problem.  Slight  modifications  in  the  initial  guess  can 
sometime  remedy  this  situation. 
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(3)  POST’s  input  can  include  as  many  as  25  control  variables  for  a  given  problem 
and  for  this  analysis  there  could  be  up  to  15.  However,  POST’s  algorithm 
becomes  even  more  sensitive  to  the  initial  guess  the  more  control  parameters 
it  is  given  to  optimize.  It  was  observed  that  the  algorithm  became  even  more 
local  as  the  number  of  control  parameters  was  increased.  For  this  reason,  it 
was  necessary  to  perform  the  optimization  in  steps.  This  was  done  by 
optimizing  the  first  few  phases  while  fixing  the  others  then  fixing  the  initial 
phases  to  those  results  and  moving  on  to  the  next  few  phases.  When  the  entire 
trajectory  had  been  through  this  process  a  few  times  the  changes  in  the  overall 
trajectory  became  small.  This  procedure,  although  necessary  for  POST  to  be 
more  global,  does  make  the  problem  more  user  intensive. 

(4)  Another  area  of  user  intervention  is  in  characterizing  the  constraints.  For 
instance,  the  heating  rate  can  be  limited  as  an  inequality  constraint  while  only 
allowable  angles  of  attack  are  used  or  the  heating  rate  can  be  maintained 
constant  using  higher  than  allowable  angles  of  attack  with  an  inequality 
constraint  on  a. 

These  are  only  some  of  the  ways  a  user  can  have  a  great  affect  on  the  outcome  of 
a  direct  method  optimization,  and  even  after  all  of  this  effort,  the  solution  is  not  optimal 
unless  it  can  be  shown  to  satisfy  those  conditions  of  the  RMP.  Thus,  an  in  depth 
knowledge  of  those  conditions  can  save  a  user  a  great  deal  of  work  by  narrowing  the 
user  s  effort.  An  example  will  be  shown  later  of  an  instance  where  the  user’s 
characterization  of  the  problem  led  to  not  obtaining  the  best  possible  solution  to  the 
problem. 

D.  TRAJECTORY  RESULTS 

The  results  obtained  in  the  investigation  of  the  two  trajectory  characterizations 
above  are  shown  below.  These  plots  represent  the  final  product  of  over  100  runs  of  the 
POST  program.  As  previously  noted  the  Aerocruise  heating  rate  profile  is  not  exactly 
constant  on  the  heat  rate  control  arc.  The  use  of  a  simulation/optimization  code  which 
allowed  the  controls  to  be  input  by  equations  vice  look  up  tables  would  remedy  this 
problem. 
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(1)  Aerobang  plots 
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Figure  17:  Velocity  vs.  Time 
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The  above  plots  show  the  entire  trajectory.  The  atmosphere  is  only  active  below 
1 1 1  km  which  corresponds  to  1271  to  3641  sec.  The  remaining  plots  will  show  this 
atmospheric  pass  only. 
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Figure  21 :  Altitude  vs.  Time 


44 


1000  1500  2000  2500  3000  3500  4000 
time  (sec) 

Figure  23:  Angle  of  Attack  vs.  Time 
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Figure  26:  Mach  Number  vs.  Time 

(2)  Aerocruise  plots 
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Figure  27:  Altitude  vs.  Time 
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Figure  28:  Velocity  vs.  Time 
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Figure  29:  Flight  Path  Angle  vs.  Time 
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Figure  30:  Propellant  Weight  vs.  Time 


The  above  plots  show  the  entire  trajectory.  Again,  the  atmosphere  is  only  active 
below  1 1 1  km  which  corresponds,  in  this  case,  to  1271  to  3641  sec.  The  remaining  plots 
will  show  this  atmospheric  pass  only. 
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Figure  36:  Thrust  vs.  Time 
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Figure  37:  Mach  Number  vs.  Time 
(3)  Comparison  of  Results 

In  [Ref.  3],  two  trajectory  characterizations  were  explored.  The  first  was  the 
previously  discussed  Aerocruise,  and  the  second  was  called  Aeroglide  by  the 
author,  but  is  not  the  same  as  the  Aeroglide  previously  discussed  in  this  thesis  as 
Aeroglide  refers  to  an  unconstrained  trajectory  in  this  thesis.  The  constrained  arc 
in  Aerocruise  of  [Ref.  3]  is  not  among  those  allowed  by  the  RMP  and  (Ref.  3j’s 
Aeroglide  trajectory  does  not  even  contain  a  constrained  arc.  This  is  an  example 
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where  the  analysis  is  incomplete  since  neither  of  the  trajectory  characterizations 
used  satisfy  the  necessary  conditions  of  the  RMP. 

When  POST  was  used  to  optimize  the  Aerocruise  trajectory,  using  the 
results  of  [Ref.  3]  as  a  starting  point,  POST  produced  a  trajectoiy  which  achieved 
24.886  degrees  of  inclination  change  compared  to  23.02  degrees  from  [Ref.  3], 
The  dynamics  and  parameterization  of  the  problem  and  vehicle  are  identical 
between  the  two  user’s  with  one  small  exception.  In  this  thesis,  the  skin  friction 
term  of  the  axial  force  coefficient  is  assumed  constant  for  a  representative 
altitude.  This  term  would  lead  to  higher  drag  coefficients  at  higher  altitudes  and 
could  account  for  some  of  the  performance  increase.  The  only  other  difference  in 
the  analyses  is  in  the  code  used.  A  non-linear  programming  code  named  GRG2  is 
used  in  [Ref.  3], 

When  POST  was  used  to  optimize  the  Aerobang  trajectory,  POST 
produced  a  trajectory  which  achieved  25.229  degrees  of  inclination  change.  This 
is  only  a  .343  degree  inclination  difference  but  does  represent  at  least  equivalent 
performance  if  not  superior  performance  to  the  Aerocruise  maneuver.  Again,  it  is 
not  implied  here  that  these  results  should  prove  the  superiority  of  either 
maneuver.  Only  that  a  performance  increase  was  obtained  when  using  a 
characterization  which  was  consistent  with  the  necessary  conditions  of  the  RMP. 
This  also  shows  that  by  using  the  knowledge  gained  from  a  theoretical 
understanding  of  the  problem,  one  can  influence  the  results  in  a  positive  way. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

From  what  has  been  discussed  and  demonstrated  it  is  clear  that  whenever  a  user 
embarks  on  using  a  direct  method  of  optimization  to  solve  a  problem,  especially  one  as 
complex  as  the  optimal  aeroassisted  orbital  plane  change,  he  must  first  ensure  that  he  has 
a  thorough  understanding  of  the  theoretical  background  of  the  problem.  Without  this 
knowledge,  the  results  of  any  analysis  will  always  be  in  question.  These  types  of 
methods  are  very  sensitive  to  the  input  given  by  the  user  and,  although  results,  will  be 
obtained,  the  validity  and  quality  of  those  results  are  most  affected  by  the  quality  of  the 
user  input. 

It  should  also  be  pointed  out  that  even  though  this  thesis  could  be  used  to  support 
the  opinion  that  some  maneuver  is  superior  to  another,  nothing  in  this  thesis  has  that 
strength.  Even  if  a  claim  of  superiority  could  be  made,  no  claim  of  optimality  could  be 
made  since  no  attempt  was  made  to  check  that  all  of  the  necessary  conditions  were  met. 
This  is  an  area  of  further  study.  Also,  even  if  optimality  had  been  shown  here,  it  would 
only  be  for  the  MRRV  from  a  very  specific  orbit  into  another  very  specific  orbit,  a  small 
subset  of  the  overall  problem  of  designing  an  optimal  vehicle  to  carry  out  missions  which 
would  utilize  this  type  of  maneuver. 

The  theory  and  numerical  analysis  performed  in  this  thesis  represent  a  good 
portion  of  what  needs  to  be  understood  when  analyzing  this  type  of  problem.  It  should 
also  be  noted  that  the  theory  is  applicable  to  any  atmospheric  trajectory  optimization  with 
a  heating  rate  constraint.  It  also  offers  some  questions  like  “How  can  two  vastly  different 
trajectory  characterizations  have  so  close  a  final  result?”,  “Does  this  mean  that  the 
optimal  change  is  close  to  that  value?”,  and  “What  would  change  in  the  results  if  the 
heating  rate  limit  were  increased  or  decreased?” 

It  is  clear  from  the  results  of  the  numerical  analysis  that  there  is  great  potential  in 
this  type  of  maneuver  for  fuel  and  mass  savings  as  well  as  pure  mission  flexibility  in  the 
design  of  future  spacecraft.  The  best  trajectory  found  in  this  thesis  results  in  a  32%  fuel 
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savings  over  the  purely  exoatmospheric  maneuver.  This  problem  offers  such  benefits 
that  it  must  continue  to  be  studied. 

B.  RECOMMENDATIONS  FOR  FURTHER  STUDY 

The  direct  extension  of  this  work  which  is  recommended  is  to  use  the  direct 
method  to  generate  a  “guess”  to  be  used  in  the  solution  to  the  TPBVP.  That  algorithm 
could  than  be  incorporated  into  a  parameter  optimization  problem  to  solve  definitively 
this  optimal  control  problem.  The  initial  guess  needed  for  the  TPBVP  solution  is  in  the 
form  of  the  initial  value  for  each  of  the  costates.  Seywald  discusses  a  method  of  costate 
estimation  in  [Ref.  15]  which  shows  great  promise  in  this  area.  This  is  also  an  area, 
however,  where  in  depth  knowledge  of  the  theoretical  underpinnings  of  the  problem  is  a 
must.  Note  that  the  discussion  in  [Ref.  15]  uses  a  Minimum  version  of  the  RMP. 

After  this  problem  is  solved  the  investigation  would  be  extended  to  other 
missions.  For  elliptical  orbits,  the  question  comes  up  as  to  where  in  the  orbit  is  best  to 
initiate  such  a  maneuver.  Then  the  problem  could  expand  into  vehicle  design  and 
mission  design.  There  is  a  vast  amount  to  be  learned  in  this  area  of  research. 
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APPENDIX  A.  SAMPLE  POST  INPUT  FILE  (AEROBANG) 


The  following  file  is  an  example  of  the  input  that  the  user  provides  to  the  POST  program 
for  optimizing  an  Aerobang  trajectory.  This  particular  file  was  used  to  optimize  the  heat 
rate  control  arcs  with  the  rest  of  the  arcs  fixed  from  previous  analysis.  Tabs  were 
inserted  into  this  file  for  aligning  the  comments,  however,  POST  will  read  tabs  as  errors. 
A  sample  MATLAB  ®  script  is  also  included  to  show  a  method  of  generating  the  tables 
required  in  the  input  format  for  POST. 


ISsearch 
srchm  =  5, 
maxitr  =  30, 
ipro  =  0, 
ioflag  =  3, 
opt  =  1.0, 
optvar  =  3hinc  , 
optph  =  100.0, 
wopt  =  0.03333, 
coneps  =  89.9, 
c 

c  Control  Variables 
c 


/accelerated  projected  gradient 
/max  number  of  iterations 

/metric  units  in-out 
/+1  for  maximization 
/optimization  variable,  inclination 
/optimization  phase 
/weighting  factor  (1/  opt  var) 
/gradient  angle  tolerance 


nindv  —  8,  /number  of  control  variables 

indvr  =  6hbnkpc2, 6hbnkpcl,  5hcritr, 
indph  =  45,  45,  50, 

u  =  -.0122452,  58.7064,  401.461, 
pert  =  .0001,  .0001,  .0001, 

c 


indvr(4)=  6hbnkpc2,  6hbnkpcl,  5hcritr, 
indph(4)  =  50,  50,  60, 

u(4)  =  .03000305,  61.4429,  440.0094, 

pert(4)  =  .0001,  .0001,  .0001, 

c 


indvr(7)—  6hbnkpc2,  6hbnkpcl, 
indph(7)  =  55,  55, 

u(7)  =  .0299218,  65.1851, 

pert(7)  =  .0001,  .0001, 

c  Target  Variables 


c 


ndepv  —  1,  /number  of  target  variables 

c 
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c  Equality  Constraints 
c 

depvr  =  6hwprop , 
depval  =  0, 

deptl  =  50, 

depph  =  100, 

c 
$ 

l$gendat 

maxtim  =200000.0,  /200000  seconds 

altmin  =  40000.0,  /40  km 

pmc  =  20 ,  /print  interval  20  sec 

pmca  =  20 , 

pine  =  20 , 

pmt(91)=  4halta,  4haltp,  3hrnc,  6henergy,  5hvcirc,  4hmass, 
pmt(97)  =  6hgenvl ,  6hgenv7 , 6hgenv3 , 6hrgenv , 
pmt(101)=  6hxmaxl ,  3hdt ,  6hgenv9 , 5hpstop, 
c 

event  =  10,  /INITIALIZE  185  KM  ORBIT 

fesn  =  100,  /final  event 

title  =  Oh*  Aerobang  with  25k  N  of  Fuel  *, 

npc(l)  =  1,  /conic  calcs 

npc(2)  =  3,  /Laplace  conic  integration  method 

npc(3)  =  5,  /initial  position 

altp  =  185.0,  /perigee  altitude 

alta  =  185.0,  /apogee  altitude 

inc  =  0.0,  /inclination 

npc(5)  =  0,  /no  atmosphere 

npc(8)  =  0,  /aerodynamic  coefF=  none 

sref  =  11.7,  /aero  reference  area 

lref  =  7.7,  /aero  reference  length 

npc(9)  =  0,  /propulsion  type,  no  thrust 

npc(  15)  =  1,  /aeroheating  rate  calculations 

npc(16)  =  1,  /spherical  planet 

omega  =  0,  /rotation  rate  of  planet 

npc(30)  =  0,  /vehicle  weight  model  to  be  used 

wgtsg  =  48041,  /vehicle  wt  (N) 

wpropi  =  25643,  /initial  propellant  wt  (N) 

iguid  =3,0,1,  /inertial  aero  angles  guidance 

monx  =  6hheatrt,  /monitor  variables 

betarg  =  4hone ,  /used  to  keep  POST  from  calculating  beta 

endphs  =  1, 

$ 

l$gendat 

event  =  20,  /DEORBIT  IMPULSE 
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/start  this  event  at  t=5  sec 


critr  =  4htime , 
value  =  5.0, 

npc(2)  =  1,  /RK-4  integrator 

npc(9)  =  3,  /Delta  V  insertion 

isdvin  =  1,  /Delta  V  input  flag 

dvimag  =  -47.07,  /magnitude  of  Delta  V 

ispv  =  300.0,  /isp  value 

endphs  =  1, 

$ 

lSgendat 

event  =  40,  /TURN  ON  ATMOSPHERE  AND  DESCEND 

critr  =  6haltito,  /when  altitude  =  120km 

value  =111000.0, 

npc(2)  =  1,  /RK-4  integrator 

npc(5)  =  2,  /1962  U.S.  standard  atmosphere 

npc(8)  =  2,  /aerocoeff,  lift  +  drag  in  tables 

npc(9)  =  0,  /no  thrust 

npc(15)  =  1,  /aeroheating  on 

iguid  =  0,1,0,  /aero  angles  guidance,  separate  channels 

iguid(6)  =  1,0,1,  /bank  and  aoa  polynomials 

alppc  =  37.50,  .002777, 

bnkpc  =  63. 11,  .0100016, 

$ 

lStblmlt 

$ 

l$tab 

/table  of  drag  coefficient  vs.  alpha 

TABLE  =6HCDT  ,  1,6HALPHA,  45,1,1,1, 

0,0.0370,  2,0.0330,  4,0.0310,  6,0.0313,  8,0.0341, 

10,0.0397, 12,0.0482, 14,0.0600,  16,0.0753, 18,0.0942, 

20,0.1171, 22,0.1442, 24,0.1757, 26,0.2119, 28,0.2530, 

30,0.2993, 32,0.3509, 34,0.4081, 36,0.4712, 38,0.5404, 

40,0.6160, 42,0.6981, 44,0.7871,  46,0.8831, 48,0.9864 
50,1.0972, 52,1.2158,  54,1.3425,  56,1.4774, 58,1.6207, 

60,1.7729, 62,1.9340, 64,2.1043, 66,2.2840, 68,2.4735 
70,2.6728, 72,2.8824, 74,3.1023, 76,3.3329, 78,3.5744 
80,3.8270,  82,4.0910,  84,4.3666,  86,4.6540,  88,4.9536, 
ixtrp  =  0,0, 0,0, 

$ 

l$tab 

/table  of  lift  coefficient  vs.  alpha 
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TABLE  =  6HCLT  ,  1, 6HALPHA ,  45,1,1,1, 

0,-0.0334,  2,-0.0211,  4,-0.0046,  6,0.0158,  8,0.0400, 

10,0.0674, 12,0.0980, 14,0.1314, 16,0.1674, 18,0.2055, 

20,0.2457, 22,0.2876, 24,0.3309, 26,0.3753, 28,0.4206, 

30,0.4664, 32,0.5126, 34,0.5587, 36,0.6046, 38,0.6500, 

40,0.6945, 42,0.7380, 44,0.7800, 46,0.8204, 48,0.8588, 

50,0.8951,  52,0.9288, 54,0.9597,  56,0.9876,  58,1.0121, 

60,1.0330, 62,1.0500, 64,1.0628, 66,1.0711, 68,1.0747, 

70,1.0733, 72,1.0666, 74,1.0543, 76,1.0361,  78,1.0118, 

80,0.9811, 82,0.9436,  84,0.8992,  86,0.8475,  88,0.7883, 
ixtrp  =  0, 0, 0, 0, 

endphs  =  1, 

$ 

lSgendat 

event  =  45,  /BEGIN  ZERO  THRUST  HEAT  CONTROL  ARC 

critr  =  6hheatrt, 

value  =11.3e+05,  /100  Btu/fitA2*s 

npc(8)  =  2,  /aero  coeff  fin  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  0,  /turn  off  rocket  engine 

npc(15)  =  1,  /aeroheating  calcs 

iguid  =  0, 1, 0,  /aero  angles  guidance 

iguid(6)=  2,0,1,  /aoa  angle  from  tables 

bnkpc  =  58.7064,  -.0122452, 

$ 

lStblmlt 

/list  of  table  multipliers 

genv7m(l)=  5.85, 6hdens  ,  /  Aref72*density 

genv6m(l)=  .15873  ,6hgenv9,  /  .5*beta/3.15 

genv4m(2)=  4hmass  ,  /  mass 

genv3m(2)=  5hgenv4 ,  /  mass*sin(gamma) 

$ 

l$tab 

table  =  5htvclt ,  0, 14679.0, 

$ 

l$tab 

table  =  6halphat,  2, 5hgenv3 , 5hgenv7 , 81, 96, 1,1, 1,1, 1,1, 1,1, 

This  was  a  bivariant  table  which  took  genv3  and  genv7  and  output  the  required  alpha. 
This  is  the  solution  to  p=0  for  the  case  of  T=0.  See  the  MATLAB  ®  script  for  now  to 
generate  this  type  of  table. 
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ixtrp  =  0,0,0,0, 

$ 

l$tab 

/table  for  vA2  multiplied  above  by  % 
table  =  6hgenv3t,  1, 5hgenv5, 2, 1, 1, 1, 
-10000.0,-10000.0,10000.0,10000.0, 

$ 

l$tab 


/table  of  sin(gamma)  multiplied  above  by  mass 
table  =  6hgenv4t,  1, 6hgammaa,  13, 1, 1, 1, 

-10.0000,  -0.17365,  -8.0000,  -0.13917,  -6.0000,  -0.10453, 

^.0000,  -0.06976,  -2.0000,  -0.03490,  -0.5000,  -0.0087, 

0.0000,  0.00000,  0.5000,  0.0087,  2.0000,  0.03490, 

4.0000,  0.06976,  6.0000,  0.10453,  8.0000,  0.13917, 

10.0000,  0.17365, 

$ 

l$tab 


/table  of  atmospheric  |3  taken  from  1962  data.  This  was 
/done  by  simulating  a  trajectoiy  the  crashed  to  the  ground 
/using  MATLAB  ®  to  convert  altitude  and  density  to  p. 
table  =  6hgenv9t,  1, 6haltito,  58, 1,  -1, 1, 

110120,  1.55e-04, 108333,  1.51e-04, 106523,  1.56e-04, 104715, 1.63e-04, 

102910,  1.69e-04, 101111, 1.76e-04,  99319, 1.74e-04,  97535,  1.78e-04, 

95762,  1.83e-04,  94002, 1.88e-04,  92257, 1.93e-04,  90531, 1.95e-04, 

88828,  1.84e-04,  87154, 1.84e-04,  85515,  1.84e-04,  83919, 1.84e-04, 

82374, 1.84e-04,  80893,  1.84e-04,  79489, 1.64e-04,  78175, 1.57e-04, 

76965,  1.53e-04,  75870, 1.50e-04,  74901, 1.47e-04,  74065, 1.45e-04, 

73366,  1.43e-04,  72804, 1.41e-04,  72374, 1.40e-04,  72067, 1.39e-04, 

71869, 1.39e-04,  71764, 1.39e-04,  71731, 1.39e-04,  71750, 1.39e-04, 

71797, 1.39e-04,  71852, 1.39e-04,  71892, 1.39e-04,  71898, 1.39e-04, 

71851,  1.39e-04,  71736, 1.39e-04,  71537, 1.38e-04,  71242, 1.37e-04, 

70952,  1.37e-04,  70661, 1.36e-04,  70359, 1.35e-04,  70119, 1.35e-04, 

69891, 1.34e-04,  69650, 1.34e-04,  69384, 1.33e-04,  69089, 1.32e-04, 

68758,  1.31e-04,  68440, 1.31e-04,  68132, 1.30e-04,  67842, 1.29e-04, 

67623, 1.29e-04,  67459, 1.29e-04,  67336,  1.28e-04,  67250, 1.28e-04, 

67201, 1.28e-04,  67190, 1.28e-04, 

$ 

l$tab 

/table  of  mu/(Rearth+altito)A2+genv6 
table  =  6hgenv5t,  2, 6hgenv6 , 6haltito,  12, 7, 
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1,1,1, 1,1, 1,1,1, 

50000, 

400,  409.65,  600,  609.65,  800,  809.65,  1000, 1009.65, 
1200,  1209.65,  1400,  1409.65,  1600,1609.65,  1800,1809.65, 
2000, 2009.65,  2200, 2209.65,  2400, 2409.65,  2600, 2609.65, 
60000, 

400,  409.62,  600,  609.62,  800,  809.62,  1000, 1009.62, 

1200. 1209.62,  1400, 1409.62,  1600, 1609.62,  1800, 1809.62, 

2000.2009.62,  2200,2209.62,  2400,2409.62,  2600,2609.62, 
70000, 

400,  409.59,  600,  609.59,  800,  809.59,  1000, 1009.59, 

1200. 1209.59,  1400, 1409.59,  1600, 1609.59,  1800, 1809.59, 

2000.2009.59,  2200,2209.59,  2400,2409.59,  2600,2609.59, 
80000, 

400,  409.56,  600,  609.56,  800,  809.56,  1000, 1009.56, 

1200. 1209.56,  1400, 1409.56,  1600, 1609.56,  1800, 1809.56, 

2000. 2009.56,  2200, 2209.56,  2400, 2409.56,  2600, 2609.56, 
90000, 

400,  409.53,  600,  609.53,  800,  809.53,  1000, 1009.53, 
1200,  1209.53,  1400,  1409.53,  1600, 1609.53,  1800, 1809.53, 
2000, 2009.53,  2200, 2209.53,  2400, 2409.53,  2600, 2609.53, 
100000, 

400,  409.50,  600,  609.50,  800,  809.50,  1000, 1009.50, 

1200.1209.50,  1400,1409.50,  1600,1609.50,  1800,1809.50, 

2000.2009.50,  2200,2209.50,  2400,2409.50,  2600,2609.50, 

110000, 

400,  409.47,  600,  609.47,  800,  809.47,  1000, 1009.47, 

1200. 1209.47,  1400, 1409.47,  1600, 1609.47,  1800, 1809.47, 

2000. 2009.47,  2200, 2209.47,  2400, 2409.47,  2600, 2609.47, 
$ 

l$tab 


/table  of  vA2  multiplied  above  by  £, 
table  =  6hgenv6t,  1, 4hvela,  6, 1, 1, 1, 

5000.0000,  25000000.000, 6000.0000,  36000000.000, 
6500.0000, 42250000.000, 7000.0000, 49000000.000, 
7500.0000,  56250000.000,  8000.0000, 64000000.000 
$ 

l$tab 


/table  of  vA2  multiplied  above  by  1/2*rho*Aref 
table  =  6hgenv7t,  1, 4hvela,  6, 1, 1, 1, 

5000.0000, 25000000.000,  6000.0000, 36000000.000, 

6500.0000, 42250000.000, 7000.0000, 49000000.000, 

7500.0000, 56250000.000,  8000.0000, 64000000.000, 
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endphs  =  1, 

$ 

l$gendat 

event  50,  /BEGIN  MAX  AOA  TRANSITION  ARC 

critr  =  6htdurp ,  /begin  after  a  certain  time  period, 

value  =  401.461, 

npc(8)  =  2,  /aero  coeff  fin  tables 

npc(9)  =  1,  /rocket  engine 

npc(22)  =  3,  /throttling  paramter  from  table 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  1,  /turn  on  rocket  engine 

npc(15)  =  1,  /aeroheating  calcs 

iguid  =  0, 1, 0,  /aero  angles  guidance 

iguid(6)  =  1, 0, 1, 

alppc  =  40,0,0,0,  /  max  aoa 

bnkpc  =  61.4429,  .03000305, 

$ 

l$tblmlt 

genv7m(l)=  5.85, 6hdens  , 
genv6m(l)=  .15873  ,6hgenv9, 
genv4m(2)=:  4hmass  , 
genv3m(2)=  5hgenv4 , 

etam  =  8.893026e-5,  /1/Tmax  to  convert  required  thrust  to  throttling 

/parameter 

$ 

I$tab 

table  =  5htvclt ,  0, 14679.0, 

$ 

l$tab 

/table  of  thrust=drag+genv3 

table  =  6hetat  ,  2, 6hdrag  ,  6hgenv3 , 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 

-100000, 

-100000,  -200000,  100000,  0, 

100000, 

-100000,  0,  100000,  200000, 

ixtrp  =  0,0,0,0, 

$ 

l$tab 

/table  of  genv5  multiplied  above  by  genv4 
table  =  6hgenv3t,  1, 5hgenv5,2, 1, 1, 1, 

-10000.0,-10000.0,10000.0,10000.0, 

$ 
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l$tab 


/table  of  sin(gaimna)  multiplied  above  by  mass 
table  =  6hgenv4t,  1, 6hgammaa,  13, 1, 1, 1, 

-10.0000,  -0.17365,  -8.0000,  -0.13917,  -6.0000,  -0.10453, 

-4.0000,  -0.06976,  -2.0000,  -0.03490,  -0.5000,  -0.0087, 

0.0000,  0.00000,  0.5000,  0.0087,  2.0000,  0.03490, 

4.0000,  0.06976,  6.0000,  0.10453,  8.0000,  0.13917, 

10.0000,  0.17365, 

$ 

l$tab 


/  see  genv9  table  previous  phase 
table  =  6hgenv9t,  1, 6haltito,  58, 1,  -1, 1, 

110120,  1.55e-04,  108333,  1.51e-04,  106523, 1.56e-04, 104715, 1.63e-04, 
102910, 1.69e-04, 101111, 1.76e-04,  99319, 1.74e-04,  97535, 1.78e-04, 
95762, 1.83e-04,  94002, 1.88e-04,  92257, 1.93e-04,  90531,  1.95e-04, 
88828, 1.84e-04,  87154, 1.84e-04,  85515,  1.84e-04,  83919, 1.84e-04, 
82374, 1.84e-04,  80893, 1.84e-04,  79489, 1.64e-04,  78175  1.57e-04, 
76965, 1.53e-04,  75870,  1.50e-04,  74901, 1.47e-04,  74065,  1.45e-04,’ 
73366,  1.43e-04,  72804, 1.41e-04,  72374, 1.40e-04,  72067, 1.39e-04, 
71869,  1.39e-04,  71764, 1.39e-04,  71731, 1.39e-04,  71750, 1.39e-04, 
71797, 1.39e-04,  71852, 1.39e-04,  71892, 1.39e-04,  71898, 1.39e-04, 
71851,  1.39e-04,  71736, 1.39e-04,  71537, 1.38e-04,  71242, 1.37e-04, 
70952,  1.37e-04,  70661, 1.36e-04,  70359,  1.35e-04,  70119, 1.35e-04, 
69891,  1.34e-04,  69650, 1.34e-04,  69384,  1.33e-04,  69089, 1.32e-04, 
68758,  1.31e-04,  68440, 1.31e-04,  68132,  1.30e-04.  67842, 1.29e-04. 
67623, 1.29e-04,  67459, 1.29e-04,  67336,  1.28e-04,  67250,  1.28e-04' 
67201, 1.28e-04,  67190,  1.28e-04, 

$ 

l$tab 

/see  genv5  table  previous  phase 
table  =  6hgenv5t,  2, 6hgenv6 , 6haltito,  12, 7, 

1,1,1, 1,1, 1,1, 1, 

50000, 

400,  409.65,  600,  609.65,  800,  809.65,  1000, 1009.65, 

1200. 1209.65,  1400, 1409.65,  1600, 1609.65,  1800, 1809.65, 

2000.2009.65,  2200,2209.65,  2400,2409.65,  2600,2609.65, 

60000, 

400,  409.62,  600,  609.62,  800,  809.62,  1000, 1009.62, 

1200. 1209.62,  1400, 1409.62,  1600, 1609.62,  1800, 1809.62, 

2000.2009.62,  2200,2209.62,  2400,2409.62,  2600,2609.62, 

70000, 

400,  409.59,  600,  609.59,  800,  809.59,  1000, 1009.59, 
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1200. 1209.59,  1400, 1409.59,  1600, 1609.59,  1800,  1809.59, 

2000. 2009.59,  2200, 2209.59,  2400, 2409.59,  2600, 2609  59’ 
80000, 

400,  409.56,  600,  609.56,  800,  809.56,  1000,  1009.56, 

1200. 1209.56,  1400, 1409.56,  1600, 1609.56,  1800, 1809.56 

2000.2009.56,  2200,2209 .56,  2400,2409.56,  2600,2609.56 
90000, 

400,  409.53,  600,  609.53,  800,  809.53,  1000,  1009.53, 

1200. 1209.53,  1400, 1409.53,  1600, 1609.53,  1800, 1809.53, 

2000.2009.53,  2200,2209.53,  2400,2409.53,  2600,2609.53’ 

100000, 

400,  409.50,  600,  609.50,  800,  809.50,  1000, 1009.50, 

1200.1209.50,  1400,1409.50,  1600,1609.50,  1800,1809.50, 

2000.2009.50,  2200,2209.50,  2400,2409.50,  2600, 2609  5o’ 

1 10000, 

400,  409.47,  600,  609.47,  800,  809.47,  1000, 1009.47, 

1200. 1209.47,  1400, 1409.47,  1600, 1609.47,  1800, 1809.47, 

2000.2009.47,  2200,2209.47,  2400,2409.47,  2600,2609  47 
$ 

l$tab 


/see  previous  phase 

table  =  6hgenv6t,  l,4hvela,  6, 1, 1, 1, 

5000.0000, 25000000.000, 6000.0000, 36000000.000, 
6500.0000, 42250000.000, 7000.0000, 49000000.000’ 
7500.0000,  56250000.000,  8000.0000, 64000000.000’ 

$ 

l$tab 


/see  previous  phase 

table  =  6hgenv7t,  1, 4hvela,  6, 1, 1, 1, 

5000.0000, 25000000.000,  6000.0000, 36000000.000, 
6500.0000, 42250000.000, 7000.0000, 49000000.000,’ 
7500.0000,  56250000.000,  8000.0000, 64000000.000’ 
endphs  =1, 

$ 

l$gendat 

event  -  55,  /BEGIN  MAX  THRUST  ARC 


critr  =  6heta  , 
value  =  1, 

npc(22)  =  0, 

/begin  when  thrust  reaches  max 

npc(8)  =  2, 

/aero  coeff  fin  tables 

npc(9)  -  1, 

/rocket  engine 

iwdf  =  2, 

/calculate  using  isp  and  thrust 

ispv  =  300.0, 

/isp 
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iengmf  =  1,  /turn  on  rocket  engine 

npc(15)  =  1,  /aeroheating  calcs 

iguid  =  0,1,0,  /aero  angles  guidance 

iguid(6)=  2,0,1,  /aoa  angle  from  tables 
bnkpc  =  65.1851,  .0299218, 

$ 

lStblmlt 

genv7m(l)=  5.85, 6hdens  , 
genv6m(l)=  .15873  ,6hgenv9, 
genv4m(2)=  4hmass  , 
genv3m(2)=  5hgenv4 , 

$ 

l$tab 

table  =  5htvclt ,  0, 14679.0, 

$ 

l$tab 

table  =  6halphat,  2, 5hgenv3 , 5hgenv7 ,111, 96, 1,1, 1,1,1, 1,1,1, 

alpha  table  as  in  phase  45  but  solution  of  p=0  for  T=Tmax 

ixtrp  =  0, 0,0,0, 

endphs  =  1, 

$ 

lSgendat 

event  =  60,  /Ascend  to  top  of  atmosphere 

critr  =  6htdurp , 
value  =  440.0094, 

npc(8)  =  2,  /aero  coeff  in  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  0,  /turn  off  rocket  engine 

npc(15)  =  1,  /aeroheating 

iguid  =  0, 1, 0,  /aero  angles  guidance 

iguid(6)  =  1,0,1, 

bnkpc  =  76.28, -.010020,  / 

alppc  =  22.96, -.01019, 

bnkarg  =  6htdurp ,  /used  here  so  bank  angle  is  continuous  through  events  60. 1 

/and  60.2 

alparg  =  6htdurp ,  /same  as  bank  angle 

$ 

lStblmlt 

$ 

IStab 

table  =  5htvclt ,  0, 14679.0, 
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endphs  =  1, 

$ 

l$gendat 

event  =  60.1,  /thrust  if  needed 

critr  =  6hgammai,  /if  don’t  get  to  1 1 1  km  then  thrust  at  apogee 

value  =0, 

npc(8)  =  2,  /aero  coefF  in  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  1,  /turn  on  rocket  engine 

npc(15)  =  1,  /aeroheating 

iguid  =  0, 1, 0,  /aero  angles  guidance 

iguid(6)  =  1,0,1, 

$ 

l$tblmlt 

$ 

l$tab 

table  =  5htvclt ,  0, 14679.0, 
endphs  =  1, 

$ 

l$gendat 

event  =  60.2,  /thrust  if  needed 

critr  =  6halta  ,  /if  don't  get  to  1 1 1  km  then  thrust  at  apogee 

value  =  185.0, 

npc(8)  =  2,  /aero  coefF in  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  0,  /turn  off  rocket  engine 

npc(15)  =  1,  /aeroheating 

iguid  =  0, 1, 0,  /aero  angles  guidance 

iguid(6)  =  1,0,1, 

$ 

l$tblmlt 

$ 

l$tab 

table  =  5htvclt ,  0, 14679.0, 
endphs  =  1, 

$ 

l$gendat 

event  =  80,  /EXIT  ATMOSPHERE 

critr  =  6haltito, 
value  =111000.0, 
npc(5)  =  0, 


/no  atmosphere 


npc(8)  = 

/no  aero  coeff 

iengmf  = 

o. 

/turn  off  rocket  engine 

iguid  = 

o,  1,0, 

/  aero  angles  guidance 

iguid(6)  = 

1,1,1, 

banki  = 

o, 

alphi  = 

o. 

betai  = 

o. 

betpc  = 

0, 0, 0, 0, 

alppc  = 

o,  o,  0, 0, 

bnkpc  = 

0, 0, 0, 0, 

endphs  = 
$ 

lSgendat 

1, 

event  = 

80.1, 

/START  REORBIT  IMPULSE 

critr  =  6hgammai, 

/Start  the  reorbit  impulse  if  gamma=0 

value  = 

o. 

npc(9)  = 

1, 

/rocket  engine 

npc(5)  = 

0, 

/no  atmosphere 

npc(8)  = 

0, 

/no  aero  coeff 

iwdf  = 

2, 

/calculate  using  isp  and  thrust 

ispv  = 

300, 

iengmf  = 

$ 

IStblmlt 

$ 

IStab 

1, 

/turn  on  rocket  engine 

table  -  5htvclt ,  0, 1500000.0,  /1000  times  normal  max  thrust 

endphs  =  1, 

$ 

lSgendat 


event  = 

80.2, 

/END  REORBIT  IMPULSE 

critr  =  6halta  , 

/stop  the  impulse  when  apogee=l  85km 

value  = 

185.1, 

npc(2)  = 

3, 

/Laplace  conic  integration  method 

dt  = 

100, 

pmc  = 

200, 

/print  interval  200  sec 

pmca  = 

200, 

pine  = 

200, 

npc(5)  = 

o. 

/no  atmosphere 

npc(8)  = 

0, 

/no  aero  coeff 

iengmf  = 

o. 

/turn  off  rocket  engine 

iwdf  = 

$ 

IStblmlt 

$ 

2, 

/calculate  using  isp  and  thrust 
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l$tab 

table  =  5htvclt ,  0,  1500000.0, 
endphs  =  1, 

$ 

l$gendat 

event  =  90,  /CIRCULARIZATION  IMPULSE 

critr  =  6haltito, 

value  =  185000.0, 

npe(9)  =  3, 

ispv  =  300, 

isdvin  =  0.0, 

endphs  =  1, 

$ 

l$gendat 

event  =  100,  /END  RUN 

critr  =  6htdurp , 

value  =  5.0, 

endphs  =  1, 

endprb  =  1, 

endjob  =  1, 

$ 

The  following  is  an  example  of  a  MATLAB  ®  script  used  to  generate  some  of  the  tables 

%  compute  the  tables  needed  for  the  heat  rate  control  portion 
%  of  the  orbit.  Makes  a  function  of  required  angle  of  attack  for  the  max  thrust  arc. 

%  this  file  was  named  MRRVTtab.m  and  executed  at  the  MATLAB  ®  command  line 

a=0; 

for  gv7=5000: 1 000: 1 00000, 
a=a+l; 
b=0; 

for  gv3=-80000: 1000:30000, 

1  j 

pol=[-gv7*4.9852e-6  -14679*.48907*(pi/180)A2-gv7*2.2429e-4  - 
14679*.00305*pi/l  80+gv7*2.4732e-3  14679*  1 .0024-gv7*3.6994e-2-gv3]; 

%  the  above  polynomial  is  for  the  solution  of  p=0  for  max  thrust  as  a  function  of  genv7 
%  and  genv3. 

alphaposs=roots(pol); 

%  the  possible  alpha’s  are  the  roots  of  the  polynomial 
%  the  below  if/then  structure  finds  the  real  root. 

if  (imag(alphaposs(  1  ))=0  &  imag(alphaposs(2))=0  &  imag(alphaposs(3))=-0) 
alphal=max(alphaposs); 
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else 


forx=l:3, 

if  imag(alphaposs(x))=0 
alpha  1  =alphaposs(x); 
end 
end 
end 

if  alphalO  %  if  the  real  root  is  negative  0  is  substituted  this  reduces 

%  interpolation  errors 
alpha  1=0; 
end 

alpha(a,b)=alphal; 

end 

end 

fid— fopen(,MRRVTtableVw');  %  opens  a  file  for  the  output 

a=0; 

%  This  section  creates  the  table  using  formatted  output  to  above  file  3  columns  of  pairs 
for  gv7=5000: 1000: 100000, 
a=a+l; 
b=-2; 

fprintf(fid,'%8.0f,\n',[gv7]); 
for  gv3=-80000:3000:30000, 
b=b+3; 
for  p=0:l, 

fprintf(fid,'%9.0f,%5.2f,',[gv3+p*  1 000  alpha(a,b+p)]); 
end 

fprintf(fid,'%9. 0f,%5.2f,\n',  [gv3+2000  alpha(a,b+2)]); 
end 
end 
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APPENDIX  B.  SAMPLE  POST  INPUT  FILE  (AEROCRUISE) 


The  following  file  is  an  example  of  the  input  that  the  user  provides  to  the  POST  program 
for  optimizing  an  Aerocruise  trajectory.  This  particular  file  was  used  to  optimize  the 
heat  rate  control,  and  atmospheric  exit  arcs  with  the  rest  of  the  arcs  fixed  from  previous 
analysis.  Tabs  were  inserted  into  this  file  for  aligning  the  comments,  however,  POST 
will  read  tabs  as  errors. 


ISsearch 
c 

srchm  =  5, 

maxitr  =  30, 

ipro  =  0, 

ioflag  =  3, 

opt  =  1.0, 

optvar  =6hinc  , 
optph  =  100.0, 
wopt  =  0.04, 

coneps  =  89.9, 

c 

nindv  =  6, 

indvr  =  6halppc2, 6hbnkpc2, 
indph  =  60,  60, 

u  =  0,  -.07, 

pert  =  .0001,  .0001, 

c 

indvr(3)  =  6hcritr ,  6halppcl,  6halppcl,  6hbnkpcl, 
indph(3)  =  60,  60,  50,  60, 

u(3)  =  1115.2,  40,  29.329,  60, 

pert(3)  =  .0001,  .0001,  .0001,  .0001, 

c 

c  Equality  Constraints 
c 

ndepv  =  5,  /number  of  target  var 

depvr  =  6hwprop , 

depval  =  0.0, 

deptl  =  50, 

depph  =  100.0, 

c 

c  Inequality  Constraints 

/  these  inequality  constraints  in  POST  are  evaluated  at  the 
/termination  of  the  pervious  phase  to  the  one  listed  for  the 


/activate  optimization  feature 
/max  number  of  iterations 

/ 

/metric  units  in-out 
/+ 1  for  maximization 
/optimization  variable,  inclination 
/optimization  phase 

/weighting  factor  (1/  opt  var) 

/number  of  control  var 
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c 


/constraint.  This  poses  an  interesting  problem  in  using  them  to 
/ensure  a  limit  is  not  violated  during  the  duration  of  a  phase  and 
not 

/just  at  the  end. 


depvr(2)  =  6hheatrt ,  6hheatrt ,  6hheatrt,  6halpha , 
depval(2)=  1.135e+06,  1.13e+06,  1.13e+06,  5, 

depph(2)  =  60,  60.1,  60.2,  80, 

idepvr(2)=  1,  1,  1,  -1, 

c 
$ 

l$gendat 

maxtim  =  20000.0,  /20000  seconds 

altmin  -  40000.0,  /40  km 

pmca  =  2.0.0, 
pine  =  20.0, 
pmc  =  20, 

pmt(91)  =  4halta,  4haltp,  3hinc,  6henergy,  5hvcirc,  4hmass, 
pmt(97)  =  6hgenvl ,  6hgenv2 , 6hgenv3 , 6hrgenv , 
pmt(101)=  6hxmaxl ,  5hpstop, 


event  =10,  /INITIALIZE  185  KM  ORBIT 

fesn  =  100, 

title  =  Oh*  Aerocruise  with  25k  N  of  Fuel  *, 
c 


npc(l)  =1, 

/conic  calcs 

npc(2)  =3, 

/Laplace  conic  integration  method 

npc(3)  =5, 

/initial  position 

altp  =  185.0, 

alta  =  185.0, 

inc  =  0.0, 

npc(5)  =0, 

/no  atmosphere 

npc(8)  =0, 

/aerodynamic  coeff  =  none 

sref  =11.7, 

/aero  reference  area 

npc(9)  =0, 

/propulsion  type,  no  thrust 

npc(15)  =  1, 

/aeroheating  rate  calculations 

npc(16)  =  1, 

/spherical  planet 

omega  =  0, 

/rotation  rate  of  planet 

npc(30)  =0, 

/vehicle  weight  model  to  be  used 

wgtsg  =48041, 

/vehicle  wt  (N) 

wpropi  =25643, 

/initial  propellant  wt  (N) 

iguid  =3,0,1, 

/inertial  aero  angles  guidance 

monx  =  6hheatrt, 

/monitor  variables 

betarg  =4hone. 

endphs  =  1, 

72 


$ 

lSgendat 

event  =  20,  /DEORBIT  IMPULSE 

critr  =  4htime, 
value  =  5.0, 

npe(2)  =1,  /RK-4  integrator 

npc(9)  =  3,  /Delta  V  insertion 

isdvin  =  1,  /Delta  V  input  flag 

dvimag  =-70.2,  /magnitude  of  the  Delta  V 

ispv  =  300.0,  /isp  value 
endphs  =1, 

$ 

lSgendat 

event  =  40,  /TURN  ON  ATMOSPHERE  &  DESCEND 

critr  =  6haltito, 
value  =111000.0, 

npc(2)  =1,  /RK-4  integrator 

npc(5)  =2,  /1 962  standard  atmosphere 

npc(8)  =  2,  /aerocoeff,  lift  +  drag  in  tables 

npc(9)  =0,  /no  thrust 

npc(15)  =  1,  /aeroheating  on 

iguid  =0,1,0,  /aero  angles  guidance,  separate  channels 

iguid(6)  =  1, 0, 1,  /bank  and  aoa  polynomials 

alppc  =  40,  0, 

bnkpc  =  50,  -.02, 

$ 

lStblmlt 

$ 

l$tab 

/table  of  drag  coefficient  vs  alpha 
TABLE  =  6HCDT  ,  1, 6HALPHA ,  45,1,1,1, 

0,0.0370,  2,0.0330,  4,0.0310,  6,0.0313,  8,0.0341, 

10,0.0397, 12,0.0482, 14,0.0600, 16,0.0753, 18,0.0942, 

20,0.1171, 22,0.1442, 24,0.1757, 26,0.2119, 28,0.2530, 

30,0.2993, 32,0.3509, 34,0.4081, 36,0.4712, 38,0.5404, 

40,0.6160, 42,0.6981, 44,0.7871, 46,0.8831, 48,0.9864, 

50,1.0972, 52,1.2158, 54,1.3425, 56,1.4774, 58,1.6207, 

60,1.7729, 62,1.9340, 64,2.1043, 66,2.2840, 68,2.4735, 

70,2.6728, 72,2.8824,  74,3.1023, 76,3.3329, 78,3.5744, 

80,3.8270,  82,4.0910,  84,4.3666,  86,4.6540,  88,4.9536, 
ixtrp  =  0, 0, 0, 0, 

$ 

l$tab 

/table  of  lift  coefficient  vs  alpha 
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TABLE  =  6HCLT  ,  1, 6HALPHA ,  45,1,1,1, 
0,-0.0334,  2,-0.0211,  4,-0.0046,  6,0.0158,  8,0.0400, 
10,0.0674, 12,0.0980, 14,0.1314, 16,0.1674, 18,0.2055, 
20,0.2457, 22,0.2876, 24,0.3309, 26,0.3753, 28,0.4206, 
30,0.4664, 32,0.5126, 34,0.5587, 36,0.6046, 38,0.6500, 
40,0.6945, 42,0.7380, 44,0.7800, 46,0.8204, 48,0.8588, 
50,0.8951, 52,0.9288,  54,0.9597,  56,0.9876, 58,1.0121,’ 
60,1.0330, 62,1.0500, 64,1.0628, 66,1.0711, 68,1.0747, 
70,1.0733,  72,1.0666,  74,1.0543,  76,1.0361,  78,1.0118, 
80,0.9811, 82,0.9436,  84,0.8992,  86,0.8475, 88,0.7883, 
ixtrp  =  0, 0, 0, 0, 

endphs  =  1, 

$ 


l$gendat 
event  =  45, 

critr  =  6hheatrt, 
value  =  1.00e+06, 

npc(8)  =  2, 

npc(9)  =  1, 

iwdf  =  2, 

ispv  =  300.0, 
iengmf  =  1, 

npc(15)  =  1, 

iguid  =  0, 1, 0, 

iguid(6)  =  1,0,0, 

bnkpc(2)  =  -.80, 

alppc  =  40,0,0,0, 

endphs  =  1, 

$ 


/make  gamma  go  to  zero 

/this  value  chosen  by  trial  and  error  to  get  heat  rate  at 
/gamma=0  to  be  max  allowed 
/aero  coeff  fm  tables 
/rocket  engine 

/calculate  using  isp  and  thrust 
/isp 

/turn  on  rocket  engine 
/aeroheating  calcs 
/aero  angles  guidance 


lSgendat 
event  =  50, 
critr  =  6hgammai , 
value  =  .006, 
npc(8)  =2, 
npc(9)  =1, 
iwdf  =  2, 
ispv  =  300.0, 
iengmf  =  1, 
npc(15)  =  1, 
npc(22)  =3, 
iguid  =  0, 1,0, 
iguid(6)  =  1, 0, 2, 

$ 

IStblmlt 


/BEGIN  AEROCRUISE  ARC 

/level  flight 
/aero  coeff  fm  tables 
/rocket  engine 

/calculate  using  isp  and  thrust 
/isp 

/turn  on  rocket  engine 
/aeroheating  calcs 
/throttling  parameter  from  table 
/aero  angles  guidance 

/bank  angle  from  tables 
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genv2m(2)  =  6hgenv4  , 
genv3m(2)  =  6hthrust , 
genv4m(2)  =  4hmass, 
etam  =  6.812453e-05, 6hdrag 
$ 

l$tab 


/table  of  lifit+genv3 

table  =  6hgenvlt,  2, 6hlift  ,  6hgenv3 , 2, 2, 1, 1, 1, 1, 1, 1, 1,  l, 

-100000, 

-100000,  -200000, 100000, 0, 

100000, 

-100000, 0, 100000, 200000, 

ixtrp  =  0,0,0,0, 

$ 

l$tab 


/table  of  cos(gamma)  multiplied  above  by  genv4 

table  =  6hgenv2t,  0,1, 

$ 

l$tab 


/table  of  sin(alpha)  multiplied  above  by  thrust 
table  =  6hgenv3t,  1, 6halpha  ,44, 1,1,1, 

0,0.00000,  1,0.01745,  2,0.03490,  3,0.05234, 

4,0.06976,  5,0.08716,  6,0.10453,  7,0.12187, 

8,0.13917,  9,0.15643,  10,0.17365,  11,0.19081, 

12,0.20791,  13,0.22495,  14,0.24192,  15,0.25882, 

16,0.27564,  17,0.29237,  18,0.30902,  19,0.32557, 

20,0.34202,  21,0.35837,  22,0.37461,  23,0.39073, 

24,0.40674,  25,0.42262,  26,0.43837,  27,0.45399, 

28,0.46947,  29,0.48481,  30,0.50000,  31,0.51504, 

32,0.52992,  33,0.54464,  34,0.55919,  35,0.57358 
36,0.58779,  37,0.60182,  38,0.61566,  39,0.62932, 

40,0.64279,  41,0.65606,  42,0.66913,  43,0.68200 

$ 

l$tab 


/table  of  mu/(Rearth+altito)A2-veliA2/(Rearth+altito) 
/multiplied  above  by  mass 

table  =  6hgenv4t,  2, 6hveli  ,  6haltito ,  16, 13, 1, 1, 1, 1, 1,1,1  l 
50000, 

5000,5.75727,  5200,5.43991,  5400,5.11011,  5600,4.76787, 

5800,4.41318,  6000,4.04604,  6200,3.66646,  6400,3.27444, 

6600,2.86997,  6800,2.45305,  7000,2.02369,  7200,1.58188," 
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7400, 1.12763,  7600,0.66093, 
55000, 

5000, 5.74530,  5200,  5.42819, 
5800, 4.40226,  6000, 4.03541, 
6600, 2.86024,  6800, 2.44365, 
7400,  1.11926,  7600,0.65293, 
60000, 

5000, 5.73337,  5200,  5.41650, 
5800, 4.39136,  6000, 4.02480, 
6600, 2.85055,  6800, 2.43428, 
7400,  1.11092,  7600,0.64494, 
65000, 

5000, 5.72146,  5200,  5.40484, 
5800, 4.38050,  6000, 4.01422, 
6600, 2.84088,  6800, 2.42493, 
7400,1.10260,  7600,0.63699, 
70000, 

5000,  5.70958,  5200,  5.39321, 
5800, 4.36966,  6000, 4.00367, 
6600,2.83124,  6800,2.41561, 
7400,1.09430,  7600,0.62905, 
75000, 

5000, 5.69774,  5200,  5.38161, 
5800, 4.35886,  6000, 3.99314, 
6600, 2.82162,  6800, 2.40632, 
7400,1.08603,  7600,0.62114, 
80000, 

5000,  5.68592,  5200, 5.37004, 
5800, 4.34808,  6000, 3.98265, 
6600, 2.81203,  6800, 2.39705, 
7400,  1.07779,  7600,0.61326, 
85000, 

5000,5.67414,  5200,5.35850, 
5800, 4.33733,  6000, 3.97218, 
6600,2.80247,  6800,2.38781, 
7400,  1.06957,  7600,0.60540, 
90000, 

5000, 5.66238,  5200, 5.34699, 
5800, 4.32660,  6000, 3.96174, 
6600, 2.79293,  6800, 2.37859, 
7400, 1.06137.,  7600,  0.59756, 
95000, 

5000,  5.65065,  5200,  5.33550, 
5800,4.31591,  6000,3.95132, 
6600, 2.78342,  6800, 2.36940, 


7800,0.18179,  8000,-0.30980, 

5400,  5.09865,  5600, 4.75667, 
6200, 3.65612,  6400, 3.26440, 
7000,2.01462,  7200,1.57316, 
7800,0.17415,  8000,-0.31705, 

5400,  5.08722,  5600, 4.74550, 
6200,3.64581,  6400,3.25439, 
7000,2.00558,  7200,1.56446, 
7800,0.16655,  8000,-0.32428, 

5400, 5.07581,  5600, 4.73436, 
6200, 3.63552,  6400, 3.24441, 
7000,  1.99657,  7200, 1.55579, 
7800,0.15896,  8000,-0.33148, 

5400,  5.06444,  5600, 4.72325, 
6200, 3.62526,  6400, 3.23445, 
7000, 1.98758,  7200, 1.54715, 
7800,0.15140,  8000,-0.33867, 

5400,5.05309,  5600,4.71217, 
6200, 3.61503,  6400, 3.22453, 
7000, 1.97862,  7200, 1.53853, 
7800,0.14386,  8000,-0.34583, 

5400, 5.04177,  5600, 4.701 12, 
6200, 3.60483,  6400, 3.21462, 
7000, 1.96968,  7200, 1.52993, 
7800,0.13634,  8000,-0.35296, 

5400, 5.03049,  5600, 4.69009, 
6200, 3.59465,  6400, 3.20475, 
7000, 1.96077,  7200, 1.52136, 
7800,0.12885,  8000,-0.36008, 

5400, 5.01923,  5600, 4.67910, 
6200,3.58450,  6400,3.19490, 
7000, 1.95189,  7200,  1.51281, 
7800,0.12138,  8000,-0.36717, 

5400, 5.00800,  5600, 4.66813, 
6200,3.57438,  6400,3.18508, 
7000, 1.94303,  7200, 1.50429, 
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7400, 1.05320,  7600,0.58974,  7800,0.11393,  8000,-0.37424 

100000, 

5000,5.63895,  5200,5.32405,  5400,4.99679,  5600,4.65719, 
5800,4.30524,  6000,3.94094,  6200,3.56429,  6400,3.17529, 
6600,2.77394,  6800,2.36024,  7000,1.93419,  7200,’  1.49579^ 
7400,  1.04505,  7600,0.58195,  7800,0 .10651,  8000-0  38129" 
105000, 

5000,5.62729,  5200,5.31262,  5400,4.98562,  5600,4.64628 
5800,4.29460,  6000,3.93058,  6200,3.55422,  6400,3.16552, 
6600,2.76448,  6800,2.35110,  7000,1.92538,  7200,  1.48732," 
7400,1.03692,  7600,0.57418,  7800,0.09910,  8000,-0  38831* 
110000, 

5000,5.61565,  5200,5.30123,  5400,4.97448,  5600,4.63540 
5800,4.28399,  6000,3.92025,  6200,3.54418,  6400,3.15577’ 
6600,2.75504,  6800,2.34198,  7000,1.91659,  7200,1.47887 ’ 
7400,  1.02882,  7600, 0.56644,  7800, 0.09173,  8000,-0.39532" 
ixtrp  =  0,0,0,0, 

$ 

l$tab 


table 


/table  of  bank  angle  as  a  function  of  genv2/genvl.  The 


/is  actually  acos(rgenv) 
table  =  6hbankt ,  1, 6hrgenv ,  201, 1,1,1, 

-1.00,  180.000000,-0.99,  171.890386,-0.98,  168.521659,-0.97,  165.930132, 

-0.96,  163.739795,-0.95,  161.805128,-0.94,  160.051556,-0.93,  158  434815’ 

-0.92,  156.926082,-0.91,  155.505352,-0.90,  154.158067,-0.89,  152.873247," 

-0.88,  151.642363,-0.87,  150.458639,-0.86,  149.316583,-0.85,  148.211669" 

-0.84,  147.140120,-0.83,  146.098738,-0.82,  145.084794,-0.81,  144.09593l" 

-0.80,  143.130102,-0.79,  142.185511,-0.78,  141.260575,-0.77,  140.353889," 
-0.76,  139.464198,-0.75,  138.590378,-0.74,  137.731416,-0.73,  136.886394 
-0.72,  136.054480, -0.71,  135.234915,-0.70,  134.427004,-0.69,  133.630109," 
-0.68,  132.843643,-0.67,  132.067065,-0.66,  131.299873,-0.65,  130  541602" 
-0.64,  129.791819,-0.63,  129.050123,-0.62,  128.316134,-0.61,  127  589503" 
-0.60,  126.869898,-0.59,  126.157008,-0.58,  125.450543,-0.57,  124  750226 
-0.56,  124.055798,-0.55,  123.367013,-0.54,  122.683639,-0.53,  122.005455," 
-0.52,  121.332251,-0.51,  120.663830,-0.50,  120.000000,-0.49,  119.340582" 
-0.48,  118.685402,-0.47,  118.034297,-0.46,  117.387108,-0.45,  116.743684," 
-0.44,  116.103881,-0.43,  115.467560,-0.42,  114.834587,-0.41  114  204835 
-0.40,  113.578178,-0.39,  112.954499,-0.38,  112.333683,-0.37,  111.715617," 
-0.36,  111.100196,-0.35,  110.487315,-0.34,  109.876874,-0.33,  109.268775," 
-0.32,  108.662925,-0.31,  108.059230,-0.30,  107.457603,-0.29,  106.857956’ 
-0.28,  106.260205,-0.27,  105.664267,-0.26,  105.070062,-0.25  104477512 
-0.24,  103.886540,-0.23,  103.297072,-0.22,  102.709033, -0.2l"  102.122352," 
-0.20,  101.536959,-0.19,  100.952784,-0.18,  100.369760,-0.17,  99.787819," 
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-0.16,  99.206896,-0.15,  98.626927,-0.14,  98.047846,-0.13,  97.469592 

-0.12,  96.892103,-0.11,  96.315316,-0.10,  95.739170,-0.09,  95.163607,’ 

-0.08,  94.588566,-0.07,  94.013987,-0.06,  93.439813,-0.05,  92.865984, 

-0.04,  92.292443,-0.03,  91.719131,-0.02,  91.145992,-0.01,  90.572967, 

0.00,  90.000000,  0.01,  89.427033,  0.02,  88.854008,  0.03,  88.280869, 
0.04,  87.707557,  0.05,  87.134016,  0.06,  86.560187,  0.07,  85.986013, 

0.08,  85.411434,  0.09,  84.836393,  0.10,  84.260830,  0.11,  83.684684’ 

0.12,  83.107897,  0.13,  82.530408,  0.14,  81.952154,  0.15,  81.373073, 

0.16,  80.793104,  0.17,  80.212181,  0.18,  79.630240,  0.19,  79.047216 

0.20,  78.463041,  0.21,  77.877648,  0.22,  77.290967,  0.23,  76.702928, 

0.24,  76.113460,  0.25,  75.522488,  0.26,  74.929938,  0.27,  74.335733, 

0.28,  73.739795,  0.29,  73.142044,  0.30,  72.542397,  0.31,  71.940770’ 

0.32,  71.337075,  0.33,  70.731225,  0.34,  70.123126,  0.35,  69.512685, 

0.36,  68.899804,  0.37,  68.284383,  0.38,  67.666317,  0.39,  67.045501, 

0.40,  66.421822,  0.41,  65.795165,  0.42,  65.165413,  0.43,  64.532440, 

0.44,  63.896119,  0.45,  63.256316,  0.46,  62.612892,  0.47,  61.965703, 

0.48,  61.314598,  0.49,  60.659418,  0.50,  60.000000,  0.51,  59.336170, 

0.52,  58.667749,  0.53,  57.994545,  0.54,  57.316361,  0.55,  56.632987’ 

0.56,  55.944202,  0.57,  55.249774,  0.58,  54.549457,  0.59,  53.842992, 

0.60,  53.130102,  0.61,  52.410497,  0.62,  51.683866,  0.63,  50.949877’ 

0.64,  50.208181,  0.65,  49.458398,  0.66,  48.700127,  0.67,  47.932935’ 

0.68,  47.156357,  0.69,  46.369891,  0.70,  45.572996,  0.71,  44.765085,’ 

0.72,  43.945520,  0.73,  43.113606,  0.74,  42.268584,  0.75,  41.409622,’ 

0.76,  40.535802,  0.77,  39.646111,  0.78,  38.739425,  0.79,  37.814489, 

0.80,  36.869898,  0.81,  35.904069,  0.82,  34.915206,  0.83,  33.901262 

0.84,  32.859880,  0.85,  31.788331,  0.86,  30.683417,  0.87,  29.54136l’ 

0.88,  28.357637,  0.89,  27.126753,  0.90,  25.841933,  0.91,  24.494648, 

0.92,  23.073918,  0.93,  21.565185,  0.94,  19.948444,  0.95,  18.194872,’ 

0.96,  16.260205,  0.97,  14.069868,  0.98,  11.478341,  0.99,  8.109614 

1.00,  0.000000, 

ixtrp  =  0,0,0,0, 

$ 

l$tab 

table  =  5htvclt ,  0, 14679.0, 

$ 

l$tab 

/table  of  thrust  as  a  function  of  alpha.  Table  is 

/l/cos(alpha)inultiplied  above  by  drag, 
table  =  6hetat  ,  1, 6halpha ,  44, 1,1,1, 

0,1.00000,  1,1.00015,  2,1.00061,  3,1.00137, 

4,1.00244,  5,1.00382,  6,1.00551,  7,1.00751, 

8,1.00983,  9,1.01247,  10,1.01543,  11,1.01872, 

12,1.02234,  13,1.02630,  14,1.03061,  15,1.03528, 

16,1.04030,  17,1.04569,  18,1.05146,  19,1.05762, 
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20,1.06418,  21,1.07114,  22,1.07853,  23,1.08636, 
24,1.09464,  25,1.10338,  26,1.11260,  27,1.12233, 
28,1.13257,  29,1.14335,  30,1.15470,  31,1.16663, 
32,1.17918,  33,1.19236,  34,1.20622,  35,1.22077, 
36,1.23607,  37,1.25214,  38,1.26902,  39,1.28676, 
40,1.30541,  41,1.32501,  42,1.34563,  43,1.36733, 
ixtrp  =  0,0,0,0, 
endphs  =  1, 

$ 

l$gendat 

event  =  60,  /Ascend  to  top  of  atmosphere 

critr  =6htdurp,  /length  of  heat  control 

npc(8)  =  2,  /aero  coeff  in  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  —  0,  /turn  off  rocket  engine 

npc(15)  =  1,  /aeroheating 

igmd  =  0,1,0,  /aero  angles  guidance 

iguid(6)  =  1,0,1, 

npc(22)  =  0,  /no  more  throttling  parameter 

eta  =  1, 
bnkarg  =6htdurp, 
alparg  =6htdurp, 

$ 

l$tblmlt 

$ 

l$tab 

table  =5htvclt,0, 14679.0, 
endphs  =  1, 

$ 

l$gendat 

event  =  60.1,  /thrust  if  needed 

critr  =  6htdurp ,  /start  thrusting  after  50  seconds 

value  =  50, 

npc(8)  =  2,  /aero  coeff  in  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  1,  /turn  on  rocket  engine 

npc(15)  =  1,  /aeroheating 

iguid  =  0, 1, 0,  /aero  angles  guidance 

iguid(6)  =  1,0,1, 

$ 

l$tblmlt 
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$ 

IStab 

table  =  5htvclt ,  0, 14679.0, 
endphs  =  1, 

$ 

ISgendat 

event  =  60.2,  /stop  thrust  when  no  longer  needed 

critr  =6halta  , 
value  =  120.0, 

npc(8)  =  2,  /aero  coeff  in  tables 

npc(9)  =  1,  /rocket  engine 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300.0,  /isp 

iengmf  =  0,  /turn  off  rocket  engine 

npc(15)  =  1,  /aeroheating 

iguid  =  0,1,0,  /aero  angles  guidance 

iguid(6)  =  1,0,1, 

$ 

l$tblmlt 

$ 

l$tab 

table  =  5htvclt ,  0, 14679.0, 
endphs  =  1, 

$ 

l$gendat 

event  =  80,  /EXIT  ATMOSPHERE 

critr  =  6haltito, 
value  =111000.0, 

npc(5)  =  0,  /no  atmosphere 

npc(8)  =  0,  /no  aero  coeff 

iengmf  =  0,  /turn  off  rocket  engine 

iguid  =  0,1,0,  /  aero  angles  guidance 

iguid(6)  =  1,1,1, 

banki  =  0, 

alphi  =  0, 

betai  =  0, 

betpc  =  0, 0,  0, 0, 

alppc  =  0,0, 0,0, 

bnkpc  =  0, 0, 0, 0, 

endphs  =  1, 

$ 

ISgendat 

event  =  80.1,  /START  REORBIT  IMPULSE 

critr  =  6hgammai,  /Start  the  impulse  if  gamma=0 

value  =  0, 
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npc(9)  =  1,  /rocket  engine 

npc(5)  =  0,  /no  atmosphere 

npc(8)  =  0,  /no  aero  coefF 

iwdf  =  2,  /calculate  using  isp  and  thrust 

ispv  =  300, 

iengmf  =  1,  /turn  on  rocket  engine 

$ 

I$tblmlt 

$ 

l$tab 

table  =  5htvclt ,  0, 1500000.0,  /1000  times  normal  max  thrust 

endphs  =  1, 

$ 

l$gendat 

event  =  80.2,  /STOP  REORBIT  IMPULSE 

critr  =  6halta  ,  /stop  the  impulse  when  apogee=l 85km 

value  =  185.1, 

npc(2)  =  3,  /Laplace  conic  integration  method 

dt  =  100, 

pmc  =  200,  /print  interval  200  sec 

pmca  =  200, 

pine  =  200, 

npc(5)  =  0,  /no  atmosphere 

npc(8)  =  0,  /no  aero  coefF 

iengmF  =  0,  /turn  oFF rocket  engine 

iwdF  =  2,  /calculate  using  isp  and  thrust 

$ 

l$tblmlt 

$ 

l$tab 

table  =  5htvclt ,  0, 1 500000.0, 
endphs  =  1, 

$ 

l$gendat 

event  =  90,  /CIRCULARIZATION  IMPULSE 

critr  =  6haltito, 

value  =  185000.0, 

npc(9)  =  3, 

ispv  =  300, 

isdvin  =  0.0, 

endphs  =  1, 

$ 

l$gendat 

event  -  100,  /END  RUN 

critr  =  6htdurp , 
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value  =  5.0, 

endphs  =  1, 

endprb  =  1, 

endjob  =  1, 

$ 
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APPENDIX  C.  MAPLE  PROGRAM  OUTPUT 


The  following  file  is  an  example  of  the  output  used  from  the  MAPLE  ®  symbolic 
manipulation  program  for  generating  the  equations  for  the  indirect  analysis.  In  this  file 
some  commands  were  used  to  actually  calculate  the  equations  while  others  were  purely 
to  produce  output  to  cut  and  paste  into  the  thesis. 
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>  with(linalg): 

Warning:  new  definition  for  norm 

Warning:  new  definition  for  trace 


>  x=vector([r,v, gamma, phi, psi,m]); 


x=[r  v  y  <j)  \|/  m\ 


>  lambda=vector([lambda[r],lambda[v],lambda[gamma],lambda[phi],lambda[psi],lambda[m]]); 

A-  =  ["A,  X  X  X  X  X 
l  r  V  y  §  \\f  m_ 

>  x:=vector([r,v, gamma, phi, psi,m,lambda[r],lambda[v],lambda[gamma],lambda[phi],lambda[psi], lamb 

>  da[m]J); 

x:=\r  vytyymX  XX  X  X  l  1 

>  u=vector([T, alpha, delta]); 

u  =  [T  a  8] 

>  dx/dt=f(x,u,t); 


>  f:=vector(14); 


>  d1:=D(alpha); 


=  f(x,u,  t ) 


/:=  array(l  ..  14,  [  ]) 
dl  ~  D(a) 


>  l:=L(alpha); 

l  :=  L(a) 

>  deltaVb:=sqrt(2*mu*(1/Ra-1/Rc)/(1-(Ra/Rc)*2*cos(gamma[f])A2))-V[f]; 

I  n  ry 


deltaVb  :=J2 


[Ra  Rc_)_  y 
Ra2  eos^y^j2  ^ 


delta  Vc:=sqrt(mu/Rc)-sqrt(2*mu  *( 1/Ra-  1/Rc)/((Rc/Ra)  A2/cos(gamma[f])A2- 1)); 


deltaVc  :=  -Jl 
.  Rc  v 


1 _ 1_ 

Ra  Rc 


2 _ \2 


RaA  cos 


>  f[1]:=v*sin(gamma); 


>  dr/dt=v*sin(gamma); 


fx  '■=  v  sin(y) 


•  =  v  sin(y) 
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>  As:=(T*cos(alpha)-d1)/m; 


>  f[2]:=As-g*sin(gamma); 


dv/dt=As-g*sin(gamma); 


As Z. cos(a)  -  D(oQ 


_  T cos(a)  -D(a)  . 

/2 - - - *“(r) 


dv  T cos(a)-D(a)  .  ,  . 

- m - ^sm(r) 


>  Ar:=(l+T*sin(alpha))*cos(delta)/m; 


Ar (L(a)  +  T_ sin(a))  cos(5) 


>  f[3]:=(Ar-g*cos(gamma)+vA2*cos(gamma)/r)A r; 

(L(a)  +  T sin(a))  cos(5) 
,  _ _ m _ 

$  ~  v 

>  v*d*gamma/dt=(Ar-g*cos(gamma)+vA2*cos(gamma)/r); 


■g  cos(y)  + 


vz  cos(y) 


vdy  (L(a)  +  Tsin(a))  cos(5)  ,  ,  v2cos(y) 
_ = - - - *«**>+ — 

>  f[4]:=v*cos(gamma)*sin(psi)/r; 

.  _  v  cos(y)  sin(v|/) 


>  d*phi/dt=v*cos(gamma)*sin(psi)/r; 


d§  v  cos(y)  sin(\{/) 
dt  r 


Aw~fl+T*sin(alpha))*sin(delta)/m; 


Aw  (L(a)  +  T_ sin(a))  sin(6) 


f[5J:=(Aw/cos(gamma)-vA2/r*cos(gamma)*cos(psi)*tan(phi))Af; 


h ' 


(L(a)  +  T sin(a))  sin(5)  v2  cos(y)  cos(v|/)  tan(<t>) 
_ m  cos(y) _  r 


v*d*psi/dt=(Aw/cos(gamma)-vA2/r*cos(gamma)*cos(psi)*tan(phi)); 


>  f[6]:=-T/lsp/gO; 


v  d \|/  (L(g)  +  T sin(a))  sin(5)  v_  cos(y)  cos(\|/)  tan(<{)) 

dt  m  cos(y)  r 


>  dm/dt=-T/lsp/g[0]; 


f  := - — 

6  IspgO 
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>  g:=mu/rA2; 


c  :=  X  v  sin(y)  +  X 
r  v 


T cos(g)  -  D(g)  [i  siG(y) 
m 


+  ■ 


? 

V  rL  j 

f  A  \ 

(L(a)  +  T_ sin(g))  cos(6)  jjl  cos(y)  vl  cos(y) 

m  2  +  r 

_ _L _  j 


+ 


¥ 


(L(g)  +  T sin(g))  sin(5)  v2  cos(y)  cos(\j/)  tan(4)) 
_ m  cos(y)  r 


r 


v  cos(y)  sm(^/) 


X  T 
m 

IspgO 


-  K 


-n 


r  j 

2t - 1 

v  max 


^  '  a  ^ 

2  — — — -  1 


) 


a 


- 1 


V  max  J 

(  u 

rcos(g)  -  D(g)  -  m  sin(y)  — 

>  fors  from  1  to  6  do;f[s+6]:=0;od: 

>  fors  from  1  to  6  do;f[s+6]:—diff(c,x[s]);od: 

>  Diff(lambda[r],t)=subs({diff(lambda[r],r)=0},f[7]); 


.w 


■  +  %v' 


JJ 


r 


X  il  sin(y)  Ky 

jLx  =-2  — — - - 

dt  r  3 


p.  cos(y)  v2  cos(y) 


v  cos(y)  sin(\}/} 


X  v  cos(y)  cos(\|/)  tan((()) 


¥ 


■  +  2 


T|  m  sin(y)  jj. 


lambdaldot 


>  Diff(lambda[v],t)=subs({diff(lambda[v],v)=0},f[8]); 


f 


zrX  =  -X  sin(y) - 2  — - — 
dt  v  r  r 


\  cos(y)  \ 


(L(g)  +  T sin(g))  cos(5)  fx  cos(y)  v2cos(y) 
m  ~  2  +  r 


J 


X  cos(y)  sin(»  X  cos(y)  cos(\(/)  tan(<J>) 

-5 - +  2  — ^ — - 

r  r 

^  (L(a)  +  T sin(g))  sin(5)  vz  cos(y)  cos(\| /)  tan(4>) 

vj/V _ m  cos(y)  r 


-2r\m  sin(y)  £  v 


lambda2dot 

>  Diff(lambda[gamma],  t)=subs({diff(lambda[gamma], gamma)=0},f[9D; 
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+ 


—  X  =  -A  vcos(y)+  — 
ot  y  r 


X  ncos(y)  y  2 


2 

|J.  sin(y)  v  sin(y) 


X  v  sin(y)  sin(\|/) 
9 


(L(a)  +  rsin(a))  sin(5)  sin(y)  ^  vz  sin(y)  cos(\|/)  tan( 4> ) 
_ m  cos(y)2 _ r _ 


/  M  9 

-r\m  cos(y)  —  +  £  v 


lamhda3dot 


>  Diff(lambda[phi],t)=subs({diff(lambda[phi],phi)=0},f[10]); 


X  vcos(y)  cos(\j/)(l  +  tan(<{>)2) 
dt  %  r 


lambda4dot 


>  Diff(lambda[psi],t)=subs({diff(Iambda[psi],psi)=0},f[11]); 


lambdaSdot 


.  A,  vcos(y)  cos(\j/)  A,  v  cos(y)  sin(\|/)  tan(<|)) 

_9  ^  =  __J> _ V _ 

9/  \|/  r  r 


>  Diff(lamb da[m ], t)=subs({ diff(lambda[m],m)=0},f[  1 2 ]); 


^  X  (T cos(g)-D(g))  A^  (L(g)  +  rsin(g))  cos(5) 


lambda6dot 


A^  (L(g)  +  r sin(g))  sin(8) 
tw2  cos(y)  v 


-•nsin(y)f-^  +  £v2 


>  Diff(L,T)=diff(c,T); 


A  cos(g)  A  sin(g)cos(5)  A  sin(g)sin(8)  A 
d v  ,  Y  |  V  _ w__2 

8T  m  mv  m  cos(y)  v  Isp  gO 


-  r|  cos(g) 

>  Diff(Diff(L,  T),  T)=subs(diff(T[max],  T)=0,  diff(diff(c,  T),  T)); 


- 1  z 


T  2- - 1 

max  T 
V  max 
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f 


dr 


■L  =  - 4 


T 

v  max 


-  1 


Z  K 

- —  +  4 - 


T 

\  max 


-1 


2f 


max 


T  ^ 

2? — 1 

max 


f 


max 


>  Diff(L,delta)=diff(c, delta); 


max 


-  1 


d  T  Y 
55 

>  Diff(L,alpha)=diff(c,  alpha); 


X  (L(a)  +  71  sin(a))  sin(5)  X  (L(a)  +  T’sin(a))  cos(8) 


¥ 


m  v 


X^-7-sm(a)-|sD(a) 

a*  m 


X  n|^L(a)]+rcos(a)]Sm(6) 


%  Xffd 
Y 


m  cos(y) v 
^L(a)  +Tcos(a)  cos(5) 


m  cos(y) v 


-2 


m  v 

(n  a  > 

k  2 - -  1 

a 

^  may 


a 


wax: 


2  — — —  -  1 
a 

^  max 


-r(|  -rsin(a)-^—  D(a) 


5 


>  Diff(Diff(L,  alpha), alpha)=diff(diff(c, alpha),  alpha); 
9  X 

d2  v 


dot" 


f 

fd2  T 

f 

T  52  ^ 

\ 

-T  cos(a)  - 

-^D(a) 

X 

2^(a) 

-  T sin(a) 

V 

{5az  X 

J  Y  ^ 

'  -i-  v 

A^az  J 

J 

cos(6) 


m 


m  v 


X 


(f  o  -X 

5Z 

\ 

2  L(a) 

-  T sin(a) 

Uaa2  J 

J 

f 


m  cos(y) v 


sin(8)  K 
- -4 - 


a 


a 

V  wax: 


-  1 


2f  a  ,  ^ 
a  I  2 - -  1 


wax: 


a 


+  2 


2  — — — -  1 
a 

v _ max 


(  5 


a  a 
oa  war 


K 


V  max  ) 
z 

z 


f  oc  ^ 
2  — — — -  1 
a 


max 


a 


max 


a 

( 

\ 

2  — — — 

- 1 

2 

a 

a 

a 

2 

- 1 

V  war 

J 

max 

a 

V 

max 

j 
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-r|  -rcos(a)-  — -D(a) 

l  Ua2  JJ 

>  D( alpha)  :=Q*( CdO+2*sin(alpha)A3); 

D(a)  :=  Q  ( CdO  +  2  sin(a)3) 

>  L(alpha):=Q*(2*sin(alpha)A2*cos(alpha)); 

2 

L(a)  :=  2  Q  sin(a)  cos(a) 

These  are  the  newtonian  values  for  the  L  and  D Junctions  for  a  flat  plate. 


>  xf:=vector([r[f],v[fJ,gamma[f],phi[f],psi[f],m[fJD; 

*"[7  7  7  7  7  7] 

>  F:=vector([cos(psi[fJ}*cos(phi[f])-cos(i[fD,r[top]-r[f],v[top]-v[f],gamma[top]-gamma[f]]); 


>  J:=-m1; 


>  nu:=vector(4); 


>  for  n  from  1  to  6  do: 


cos  r  -r „  v  -v„  v  -v 
\f)  top  f  top  f  fop  j 


J  :=  -ml 


v  :=  array(  1  ..  4,  [  ]) 


>  lambdatf[n]:=-diff(J,xf[n]): 


>  for  s  from  1  to  4  do: 


>  lambdatf[n]:=lambdatf[n]+nu[s]*diff(F[s],xf[n]): 


>  for  n  from  1  to  6  do:  lambdatf[n];od; 


-V1  cos^)  sin(<y 

— v  j  sin(vy)  c°s(4y) 


Now  Iwil  try  to  look  at  the  jump  conditions  from  constrained  to  unconstrained.  The  applicable  condition  to 
transition  form  unconstrained  to  constrained  is  continuity  in  lambda. 
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>  for  n  from  1  to  6  do 


>  lamb daa[n]=lamb dab[n ]-eta  *diff(G,x[n ]);o d; 

lambdaa  ^  =  lambdab  ^  +  r\  C  pO  $  N 

lambdaa  =lambdab  _.?Cp0e(  - 

2  2  v 

lambdaa  =  lambdab 

3  3 

lambdaa  =  lambdab 

4  4 

lambdaa  =  lambdab 

5  5 

lambdaa  =  lambdab 

6  6 
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